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F.  W.  Wendt  for  his  encouragement  and  support  throughout  the  investigation. 
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ABSTRACT 


The  dynamic  response  of  a  metallic  body  when  impacted  by  a  hyper¬ 
velocity  projectile  is  a  complicated  phenomenon  involving  the  introduction 
of  damage  mechanisms  which  are  not  completely  understood,  even  when 
independently  operative.  In  Part  I  of  this  report,  tentative  phenomeno¬ 
logical  relations  are  proposed  for  the  hypervelocity  impact  regime. 

They  include  (a)  the  flow  resistance  coefficient  (which  is  of  fundamental 
importance  in  the  visco-plastic  model  for  high  speed  impact),  (b)  the 
equation  of  state  for  both  compressive  and  rarefaction  regions  of  the 
metal,  and  (c)  a  fracture  criterion  for  high  intensity  tri-axial  stress 
distributions  of  extremely  short  duration. 

In  Part  II  of  the  report  a  step-by-step  description  of  a  computational 
procedure  is  presented  for  the  solution  of  the  visco-plastic  equations 
governing  the  axisymmetric  hypervelocity  impact  situation.  The  pro¬ 
cedure  incorporates  the  above  phenomenological  relations  which  are 
necessary  for  a  realistic  description  of  the  penetration  and  cratering 
processes.  Finally,  a  brief  resume  is  given  of  the  considerations  that 
were  involved  in  developing  PIC  WICK,  an  IBM  7090  computer  program 
for  carrying  out  the  complex  and  lengthy  computational  procedure.  This 
unique  program  is  a  basic  tool  in  the  development  of  a  more  complete 
theoretical  understanding  of  the  hypervelocity  impact  phenomenon. 
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INTRODUCTION 


The  dynamic  response  of  a  metallic  body  in  impact  by  a  hypervelocity 
projectile  is  a  complicated  phenomenon  involving  the  introduction  of  damage 
mechanisms  which  are  not  completely  understood  even  when  independently 
operative.  However,  a  qualitative  model  of  crater  formation  has  evolved, 
based  on  experimental  investigations  of  impact  at  velocities  as  high  as  2 
cm/ n  sec.  *  Tests  in  which  the  actual  cratering  process  has  been  monitored 
by  high  speed  photography  and  flash  X-ray  photography  have  shown  that  the 
crater  may  continue  to  enlarge  for  several  hundred  microseconds  even  though 
only  five  to  ten  microseconds  are  required  to  use  up  the  projectile  (Reference 
1).  The  mechanism  of  crater  formation  is  therefore  essentially  one  of  cavi¬ 
tation,  the  size  and  shape  of  the  final  crater  being  determined  by  the  pressure 
wave  established  during  the  first  five  to  ten  microseconds  by  the  action  of  the 
projectile  on  the  target,  and  by  the  resistance  of  the  target  material  to  flow. 
The  flow  continues  until  the  amplitude  of  the  wave  decreases  below  the  in¬ 
trinsic  yield  strength  of  the  target  material. 

In  previous  reports  (References  2  through  4)  the  available  experimental 
data  have  been  studied  and  a  visco -plastic  model  has  been  formulated  which 
takes  into  account  the  above  characteristic  effects.  This  model  includes  the 
viscous  and  strength  effects  of  the  impacted  metals  as  well  as  the  inertial 
effect.  The  importance  of  including  all  three  effects  in  the  model  was  demon¬ 
strated  quantitatively  by  numerical  solution  of  the  system  of  equations  govern¬ 
ing  the  model  in  the  case  of  a  one -dimensional  impact  situation.  The  results 
of  these  earlier  investigations  were  discussed  in  detail  in  References  5  and  6. 

The  object  of  this  phase  of  the  theoretical  investigation  of  the  hyperve¬ 
locity  impact  phenomenon  is  the  solution  of  the  governing  equations  for  the 
axisymmetric  case  depicted  in  Figure  1.  This  configuration  requires  the 
solution  of  the  system  of  two-dimensional,  time  dependent  partial  differential 
equations  given  in  Appendix  A;  a  large  difference  exists  in  the  difficulty  be¬ 
tween  the  1-D  and  2-D  problems. 


♦The  gram -centimeter  microsecond  system  of  units  is  used  throughout  this 
report.  The  unit  of  stress  becomes  the  megabar  (mb). 
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In  developing  a  finite  difference  formulation  of  the  axiaymmetic  impact 
problem  it  ie  natural  to  attempt  an  extension  of  an  existing  scheme  devised 
for  two-dimensional  hydrodynamics;  this  has  been  our  approach  (References 
7  and  8).  In  addition  to  the  obvious  difficulties  of  extending  the  computational 
procedure  to  the  much  more  complicated  system  of  equations  governing  the 
visco -plastic  model,  the  description  of  the  axisymmetric  impact  problem  must 
also  include  provisions  for  negative  pressure  and  material  fracture  (Reference 
9).  The  fundamental  importance  of  these  phenomena  in  the  cratering  process 
is  depicted  by  the  sequence  of  illustrations  in  Figure  2.  The  rupture  of  the 
material  and  its  ejection  from  the  crater  during  the  cavitation  process  accounts 
for  a  large  percentage  of  the  final  crater  volume.  In  Part  I  of  this  report 
phenomenological  relations  are  proposed  which  permit  these  phenomena  to  be 
included  in  the  theoretical  model. 

Part  II  of  this  report  gives  a  comprehensive  summary  of  the  various  con¬ 
siderations  that  have  gone  into  the  development  of  a  numerical  scheme  for  the 
solution  of  the  axisymmetric  hypervelocity  impact  problem.  The  description 
is  complete  enough  that  a  programmer  may  work  directly  from  this  report. 
Indeed  a  computer  code  for  the  IBM  7090  (or  7094)  has  been  completed  and  is 
in  operation.  The  program  is  named  PICWICK;  its  salient  features  are 
described  in  the  final  section  of  this  report. 
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PART  I:  PHENOMENOLOGICAL  RELATIONS 


For  each  material  requiring  calculations  the  equation  of  state,  the  flow* 
resistance  coefficient,  and  the  fracture  criteria  for  a  dynamic  tri- axial  stress 
condition  must  be  specified.  Actually,  none  of  these  have  been  firmly  estab¬ 
lished  by  experiments  under  the  severe  conditions  of  pressure,  strain- rate, 
and  temperature  which  occur  during  hypervelocity  impact.  At  the  present 
state  of  knowledge  it  is  necessary  to  extrapolate  boldly  from  data  observed 
under  far  less  severe  conditions.  This  is  not  so  important  from  the  theoreti¬ 
cal  viewpoint  as  the  ability  of  the  analysis  to  include  and  compare  various 
proposed  constitutive  relations.  To  facilitate  the  evaluation  of  these  various 
hypotheses,  the  equation  of  state,  the  flow- resistance  coefficient,  and  the 
fracture  criterion  have  been  written  into  the  program  as  sub-routines  so  that 
they  may  be  readily  changed. 

In  this  part  of  the  report  tentative  constitutive  relations  are  proposed  for 
the  hype rvelocity  regime.  As  more  experimental  data  become  available,  and 
as  computational  experiments  are  performed,  improved  hypotheses  are  ex¬ 
pected  to  evolve. 


FLOW- RESISTANCE  COEFFICIENT 

The  response  of  the  projectile  and  target  materials  when  subjected  to 
hypervelocity  impact  is  governed  by  the  system  of  visco-plastic  equations 
formulated  earlier  in  this  study  (References  2  and  6).  The  equations  were 
constructed  to  bridge  the  transistion  from  the  plastic  to  the  hydrodynamic 
regimes  in  the  simplest  manner  consistent  with  the  physical  requirement 
that  the  inertial  (compressibility),  viscous,  and  strength  effects  be  included. 

Thus  the  stress  tensor  r  and  the  strain-rate  tensor,  D.^,  were  as¬ 
sumed  related  acording  to 

rik =  "  Dik  - -r M  diy "  *ik‘  “» 

the  essential  new  feature  of  the  model  was  the  introduction  of  the  flow- 
resistance  coefficient 
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M(D) 


=  oo 


.  . 

<M  a  so) 

(M  <V 


into  the  usual  momentum  and  energy  equations  of  continuum  mechanics.  The 
quantity  denotes  the  second  invariant  of  the  deformation  strain-rate  tensor, 
denotes  the  static  yield  shear  stress  of  the  material,  n  denotes  the  vis- 
cosity  factor,  and  r4  =  /i4  D4  is  a  measure  of  the  deformation  experienced  by 

the  medium.  Thus,  the  material  was  considered  rigid  if  stressed  below  its 
yield  strength,  whereas  above  this  value  it  was  assumed  to  behave  as  a 
Newtonian  viscous  liquid. 

Recently  (Reference  8),  a  small  parameter  (  >  0  was  incorporated  into 
the  definition  of  the  flow  resistance  coefficient, 

S 

M(D)  =  VQ  +  — — -  (2) 

j/d*  +  t 


It  was  introduced  chiefly  to  remove  the  moving  surface  of  separation  between 
the  rigid  and  fluid  regions  of  the  medium.  This  not  only  simplifies  the  calcu¬ 
lations,  however,  but  is  also  more  realistic  because  the  shearing  stresses 
now  depend  upon  the  strain- rate  in  a  continuous  manner.  Under  prolonged 
loading  this  model  permits  deformation  to  occur  even  for  I  r|  <  SQ,  but  the 
impact  mechanism  is  completed  long  before  such  creep  effects  can  occur. 

Actually  both  r)  0  and  SQ  are  not  constant  but  depend  on  the  thermodynamic 
state  of  the  medium: 


i?0  -**»(!.  P)  So-*S(I,  p). 

Both  decrease  in  value  if  the  specific  internal  energy  (essentially  the  temper¬ 
ature)  is  increased  while  holding  the  pressure  constant, 
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both  increase  in  value  if  the  pressure  is  increased  while  the  internal  energy 
is  held  constant, 


In  order  to  explicitly  account  for  these  effects,  the  definition  of  the  flow 
resistance  coefficients  will  be  written  as  follows:  41 

ft  s  ft  (l,  p,  D) 


2  *(I.  p)  + 


S(I,  p) 

Vd*  +  f 


(5) 


The  internal  energy,  pressure,  and  strain-rate  dependence  of  p  are  depicted 
schematically  in  Figure  3.  The  governing  equations  are  given  in  Appendix  A. 

To  appreciate  the  significance  of  the  dependence  of  ft  on  the  thermo¬ 
dynamic  variables,  substitute  the  continuity  equation  into  the  adiabatic  energy 
equation  to  obtain 

dl  -1  T  2  -1  dp  I 

#  ■  '  at" J  '•  <6> 

the  specific  internal  energy  therefore  increases  monotonically  with  increasing 
strain-rate  (compressive).  In  view  of  (3)  and  (5)  it  follows  that,  for  a  given 
final  compression,  the  flow-resistance  coefficient  of  the  medium  will  decrease 
rapidly  with  an  increase  in  the  strain-rate.  The  decrease  in  the  flow- 
resistance,  in  turn,  permits  a  further  increase  in  the  strain-rate.  The  de¬ 
creased  resistance  to  deformation  at  extreme  strain-rates  has  been  observed 
to  be  extremely  important  in  the  cratering  process,  as  it  permits  the  cavita¬ 
tion  flow  to  continue  a  long  time  after  impact  (Reference  1).  Since  the  greater 
part  of  the  crater  is  formed  during  the  cavitation  flow,  the  dependence  of  tf 
and  S  on  the  thermodynamic  state  of  the  medium  cannot  be  safely  neglected. 

This  feedback  between  the  temperature  and  further  decrease  in  flow- 
resistance  is  currently  the  subject  of  research  by  Dr.  I.  J.  Gruntfeat  (Refer¬ 
ence  10). 


The  flow-resistance  coefficient  is  not  restricted  to  the  form  in  relation  (5). 
The  numerical  method  can  treat  any  general  function  fi~  ft{ I,  p,  D). 
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SURVEY  OF  FLOW -RESISTANCE  DATA 


There  are  apparently  no  measurements  of  the  flow-resistance  coefficients 
of  materials  under  the  conditions  of  temperature,  pressure,  and  strain-rate 
that  occur  in  hypervelocity  impact.  Data  are  available,  however,  in  the 
following  isolated  thermodynamic  ranges: 

a.  Viscosity  measurements  of  molten  metals  at  atmospheric  pressure. 

b.  Longtime,  low-load  creep. 

c.  Plastic  wave  propagation  in  metal  bars. 

d.  Static  yield  stress  measurements  as  functions  of  pressure  and 
temperature. 

These  data  have  been  surveyed  to  give  rough  estimates  for  the  functions 
T)(I,  p)  and  S  (I ,  p  ). 


a)  .Molten  Metals 

Above  its  melting  temperature  the  yield  strength  of  a  metal  vanishes, 

S  =  0,  and  the  flow-resistance  coefficient  reduces  to  the  ordinary  viscosity 
factor,  n  «  T),  Experimental  values  of  p  are  therefore  available  in  the 
elevated  temperature  range  (at  atmospheric  pressure)  from  measurements 
of  the  flow-resistance  of  the  liquid  phase.  The  experimental  results  for  a 
number  of  metals  are  presented  in  Figure  4.  The  viscosity  factor  is  seen 
to  vary  with  the  absolute  temperature  according  to  the  relation 


M  o  =  A  exp  (  C/T)  . 


(7) 


The  values  of  A  and  C  are  given  in  Table  I  for  the  various  metals. 

It  has  been  observed  that  the  flow-resistance  of  a  metal  suddenly  in¬ 
creases  as  the  temperature  is  lowered  to  its  melting  point  (Reference  14). 
According  to  the  present  hypothesis,  equation  (5),  this  effect  may  be  as¬ 
sumed  to  be  a  consequence  of  S  assuming  a  non-zero  value  at  the  freezing 
temperature  rather  than  indicating  a  discontinuity  in  rf.  With  this  interpreta¬ 
tion,  the  above  exponential  relation  may  be  used  to  make  zero -order  esti¬ 
mates  for  the  solid  phase  values  of  »?. 
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The  specific  internal  energy  will  be  set  equal  to  zero  when  T  =  300°K 
and  p  =  0.  Then  equation  (7)  is  essentially  equivalent  to 


Vil,  0) 


(8) 


where  li  0  =  A,  H  \  =  cC,  H  i~  300c  and  c  is  a  mean  value  of  the  specific  heat 
at  zero  pressure.  The  values  of  not  n\,  and  2  are  in  Table  I,  as 

are  the  estimated  values  of 

T)(0,  0)  =  fi  exp(  M,  /  M  . 

o  1  c 


h)  Creep  Data 

Creep  data  is  ordinarily  generated  by  studying  the  long-term  elongation 
of  a  long  rod  subjected  to  uniaxial  tension  (along  z-axis).  Then 


r 

rr 


=  0 


P  = 


3 


r 

.  zz 


D  =  0 
rz 


div  u  s  0. 


The  flow -resistance  coefficient  and  the  axial  stress  reduce  to 


In  this  range  therefore,  H  will  be  very  large  since  <  and  dv/dz  are  small. 
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For  example,  in  the  case  of  high-purity  polycrystalline  aluminum  sub¬ 
jected  to  long  term  creep  (Reference  15): 

T  =  478°K(I  =  1.  71  x  lO-3),  r  =  2.  1  x  lO-4,  =  10"U, 

zz  o  z 

-4  * 

S  =  2  x  10  . 

Hence,  using  (10), 


-4 


2.  1  x  10 

*  ~ - rn 

3  x  10 


=  0.  7  x  10 


The  first  term  on  the  right  side  of  (9)  is  approximated  by  6  x  10"^.  It  is 
therefore  negligible  and  we  may  write  0.  7  x  10^  aj  2  x  10"^  c,  whence, 

Aluminum:  <  »  3  x  10 


Actually,  the  true  value  of  <  is  unimportant  in  impact  calculations.  It 
will  therefore  be  assigned  as  convenient  in  the  calculations  without  regard  to 
values  estimated  from  creep  data.  It  must,  however,  be  kept  small. 


c)  Plastic  Wave  Propagation 

Dynamic  stress -strain  curves  for  high -purity,  polycrystalline  aluminum 
have  recently  been  presented  by  Hauser,  Simmons,  and  Dorn  (Reference  16). 
For  static  tests  made  at  T0  =  295°K  the  yield  stress  was  VT  SQ  «  0.  77  x 
10~3.  The  dynamic  stress  at  a  total  strain  of  16%  was  observed  to  be  1.  96  x 
10"^  and  i.  36  x  10"3  for  strain  rates  of  10"^  and  10"^  respectively.  Equa¬ 
tions  (9)  and  (10)  may  again  be  used  to  estimate  t|  since  div  u  =  0  in  these 
tests: 


*The  gram-centimeter-microsecond  system  of  units  is  used  throughout  this 
report. 
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or 


Aluminum: 


r]  =  0. 04 


n  =  0. 20 

O 


for  strain -rate  of  10 


for  strain -rate  of  10 


These  results  simply  verify  that  the  viscosity  factor  cannot  be  regarded 
as  constant.  The  decrease  in  its  value  with  increase  in  the  strain -rate  may 
be  assigned  to  the  greater  specific  internal  energy  of  the  medium.  If  the 
strain -rate  10"2  is  great  enough  for  a  truly  adiabatic  situation,  the  increase 
in  specific  internal  energy  is  roughly  estimated  by 


=  (£)„ 


where  t*  =  time  duration  of  test  =  0.  16/ 10"2  =  16  ftsec.  and  (6)  may  be 
used  to  approximate  dl/dt.  This  procedure  yields  an  estimate  of  Al  sr  7.  5  x 
10’5  which  corresponds  to  a  temperature  rise  of  AT  *  8.  5°K.  If  we  assume 
that  in  this  regime  the  temperature  dependence  of  *1  is  also  of  the  form 
A  exp  (C/T)  the  estimates  for  the  parameters  would  be  A  *  0.  113  and  C  *  170 
as  opposed  to  the  values  A  =  10"*  *  and  C  =  7580  deduced  from  the  molten 
metal  data.  The  former  values  may  be  regarded  as  a  first  order  approxi¬ 
mation. 

The  values  of  VQ  listed  in  (11)  compare  favorably  with  the  value  of  0.  23 
assumed  by  Malvern  (Reference  17)  in  his  original  study  of  the  strain-rate 
effect  in  aluminum  rods.  Values  for  steel  and  copper  have  been  estimated  by 
Perzyna  (Reference  18): 


Iron: 


Copper: 


=  0.  8  gm.  cm  *  ft  sec  1  (Perzyna) 
-1  -1 

7  =  0, 4  gm.  cm  ft  sec  (Perzyna) 


Aluminum:  =  0.  23  gm.  cm  *  ft  sec  1  (Malvern). 
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~  rafa  - 


d)  Thermodynamic  Dependence  of  Yield  Strength 

The  literature  contains  sufficient  data  relating  to  the  temperature  depend¬ 
ence  of  the  yield  strength  of  metals  under  atmospheric  pressure  conditions 
(e.  g. ,  Reference  19).  The  function 


S(I.  0) 


(14) 


may  be  fitted  to  the  data  for  any  particular  metal.  No  special  significance  is 
attached  to  the  form;  it  merely  has  the  property  of  decreasing  with  an  increase 
in  1  and  contains  a  sufficient  number  of  material  parameters  to  permit  the  re¬ 
quired  flexibility. 

As  a  result  of  static  tension  tests  under  hydrostatic  pressure  Bridgeman 
concludes  (Reference  20,  p  70)  that  the  yield  stress  is  raised  by  an  amount 
roughly  proportional  to  the  pressure.  The  combined  effects  of  pressure  and 
temperature  have  apparently  not  been  investigated  even  under  static  condi¬ 
tions. 


» 


EQUATION  OF  STATE 

In  the  compressive  regime  the  most  reliable  equation  of  state  available 
has  been  determined  by  the  Los  Alamos  group  from  velocity  measurements  of 
shock  waves  induced  by  high  explosives.  The  pressure  is  expressed  as  a 
function  of  density  and  specific  internal  energy  (Reference  21). 

p  =  f(p.  I) 


where 


*2a2) 

+  *  j  [bo+(bl  +  b2  £)  <]  [  1+  *2  A(1  "  X)]  (15) 

+  ♦Cc.+  CjO  [l  -  *2<1  -M2]  j  -  U2<a2  +  b2*)J  . 
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Here,  there  are  introduced  the  notations 


<  - 


*  -  'a1 


A 


* 

0O 


(16) 


and  the  material  parameters 

a i »  b  ,  b. ,  b  ,  c  ,  c.,  <p  i 

1  2  o  1  2  o  1  *o  *1  *2 


(17) 


The  specific  internal  energy  in  the  undisturbed  material  is  zero,  1Q  =  0. 
Throughout  the  calculations  1^0. 

It  is  also  necessary  to  provide  an  equation  of  state  for  the  tensile  regime 
since  rarefaction  regions  occur  near  the  edge  of  the  projectile -target  interface 
during  the  early  stages  of  the  process,  and  near  the  lip  of  the  forming  crater 
during  the  later  stages.  Apparently  no  experimental  data  are  available  under 
these  extreme  conditions.  It  seems  reasonable  and  expedient  to  use  the 
tangent  line  of  f  ( p ,  I)  to  extend  the  equation  of  state  into  the  tensile  regime 
(Figure  5): 

mo.1)  ■  <[i?]fi0  +  M<r=o 

or 

h«’-I)  ’  -jrhr0  [»,<<* -f2*2>  +  *}[bo  +  bi<][1+  V(I-X)] 

+  0  (cQ  +  Cj  o  [l  -  f2(l  -  A)2]  J  -  ex(2C+  l)(a2+b2*)J 

(18) 

The  equation  of  state  is  therefore  assumed  to  be  given  by  the  following 
composite  formula: 


P  =  F(p,I) 
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where 


F (p,  I) 


Up,  I)  if  <  >  0 

h <p,  I)  if  C  <  0  • 


(19) 


The  material  parameters  (17)  are  listed  in  Table  II  for  aluminum  and  iron. 
Values  are  also  available  for  most  other  metals  of  interest  but  they  are 
apparently  still  classified. 


FRACTURE  CRITERION 

To  account  for  the  ejection  of  material  from  the  lip  of  the  crater  during 
formation,  it  is  necessary  to  incorporate  a  dynamic  fracture  criterion  into 
the  calculations.  If  no  such  provision  is  made,  one  would  be  tacitly  assuming 
extremely  large  fracture  strengths  for  the  projectile -target  materials.  The 
consequence  would  be  that  no  material  would  be  ejected  from  the  target  (which 
would  be  contrary  to  experimental  observations)  and  the  calculated  crater 
dimensions  would  be  much  smaller  than  are  actually  found. 


The  material  is  subjected  to  tri -axial  stress  under  extreme  conditions  of 
temperature  and  strain-rate.  No  fracture  criterion  has  been  established  under 
these  conditions.  In  the  present  formulation  it  is  assumed  that  fracture  oc¬ 
curs  whenever  both  the  pressure  is  negative  (i.  e.  ,  hydrostatic  tension  exists) 
and  the  maximum  deviatoric  stress  exceeds  a  critical  value: 


p  <  0  and 


V  £  Q 
max  cr 


(20) 


where  the  deviatoric  stress  components  are  defined  in  equation  (A-ll)  and 


V  =  max 
max 


[f 


rr 


ee 


zz 


I] 


(21) 


The  critical  value  of  acr  will  depend  upon  the  temperature  and  the  length  of 
time  during  which  the  stress  is  applied. 


A  rational  assumption  is  that  the  damage  suffered  by  a  segment  of  the 
medium  is  cumulative;  that  is,  in  each  small  time  increment  the  fracture  will 
proceed  at  a  rate  appropriate  to  the  stress  distribution  and  temperature 
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occurring  during  that  time  increment.  In  the  present  numerical  scheme, 
however,  the  stress  field  is  expressed  in  Eulerian  coordinates,  and  it  would 
be  extremely  difficult  to  account  for  the  cumulative  damage  suffered  by  the 
material  particles.  It  is  therefore  assumed  that  damage  is  cumulated 
only  for  the  time  interval,  8t,  corresponding  to  one  time  cycle  of  the  numer¬ 
ical  scheme. 


In  the  case  of  uniaxial  loading  the  fracture  criterion  most  generally  ac¬ 
cepted  for  dynamic  conditions  takes  the  form  (Reference  22): 


(22) 


where  R  is  the  gas  constant,  tj  is  the  time  in  which  the  stress  becomes 
tensile,  1 2  is  the  time  at  which  fracture  occurs,  and  tQ,  UQ,  and  y  are 
material  constants.  To  extend  this  relation  in  accordance  with  the  above 
discussion,  replace  a  by  the  critical  value  of  Vmax,  o cr,  and  cumulate 
only  for  time  8 1  to  get 


exp 


X  -  <ao 
_ cr 

T 


U> 


y_ 

R 


or,  taking  the  natural  logarithm  of  both  sides  and  rearranging  terms, 


a  * 
cr 


+[ 


X  -  T  tn(8 t/t  ) 
o 


1- 


(23) 


In  terms  of  the  specific  internal  energy,  the  relation  is  approximated  by 


a  * 

cr 


-  [x  *  &  •  ™) f*  f  ] ' 


(24) 


where  c  is  a  mean  value  of  the  specific  heat. 

From  experimental  results  over  a  wide  range  of  stresses  and  temper¬ 
atures  Zhurkov  (Reference  22,  23,  24)  has  deduced  empirical  values  for  the 
constants  tQ,  UQ,  and  y  for  a  number  of  metals.  UQ  was  found  to  be  nearly 
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equal  to  the  sublimation  energy  while  tQ  was  found  to  be  of  the  order  of  mag¬ 
nitude  of  the  vibration  period  of  the  crystal  lattice.  Zhurkov's  results  for 
platinum  and  aluminum  have  been  replotted  in  the  gram -centimeter - 
microsecond  system  of  units  and  are  given  in  Figures  6  and  7  respectively. 
The  approximate  values  of  tc,  X  and  e>  obtained  from  these  plots  are  pre¬ 
sented  in  Table  III  along  with  estimated  values  for  several  other  metals. 
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PART  II:  METHOD  OF  SOLUTION 


In  developing  a  finite  difference  formulation  of  the  equations  governing 
the  axi symmetric  impact  problem  (Appendix  A),  the  extension  of  an  existing 
scheme  devised  for  two-dimensional  hydrodynamics  is  a  natural  approach. 
Several  methods  of  treatment  have  been  used  for  those  problems  dependent 
upon  two  or  more  space  coordinates.  These  variations  usually  employ  (a) 
Lagrangian  coordinates  in  which  the  mesh  of  cells  is  imbedded  in  the  medium 
and  moves  with  it  (References  29  through  32),  (b)  Eulerian  coordinates  which 
are  not  fixed  in  the  medium  but  are  usually  stationary  in  the  laboratory  frame 
of  reference,  or  (c)  a  mixed  Euler- Lagrange  system  which  attempts  to  take 
advantage  of  the  better  features  of  both  fixed  and  moveable  coordinates  (Ref¬ 
erences  33  through  36). 

The  chief  difficulty  with  schemes  employing  Lagrangian  coordinates  is 
the  large  distortion  which  is  involved  in  the  present  problem.  The  Eulerian 
systems  have  the  disadvantage  that  to  account  for  the  projectile -target  inter¬ 
face  and  the  free  surfaces  of  the  projectile  and  target  is  extremely  difficult. 
These  are  the  basic  reasons  for  the  decision  to  adopt  the  particle -in-cell 
method  (abbreviated  PIC)  which  has  been  developed  at  Los  Alamos  by  Dr. 
Francis  H.  Harlow  and  his  colleagues. 

In  devising  the  present  formulation  the  basic  references  have  been  the 
series  of  Los  Alamos  reports  which  describe  the  PIC  method  (References  33 
through  35).  An  attempt  has  been  made  to  derive  maximum  benefit  from  the 
experiences  of  the  Los  Alamos  group  in  applying  the  method  and  to  extrapo¬ 
late  their  results  and  heuristic  concepts  to  the  more  complex  problem  at  hand. 
This  is  especially  desirable  in  using  the  PIC  method  as  it  is  still  an  experi¬ 
mental  procedure  which  does  not  have  a  firm  mathematical  basis.  The  author 
has  had  the  benefit  of  extremely  fruitful  discussions  with  Dr.  Harlow  at  Los 
Alamos. 

A  brief  description  of  the  essential  features  of  the  method  may  be  given. 
The  space  occupied  by  the  medium  is  divided  into  a  suitably  chosen  mesh  of 
fixed  cells  through  which  the  medium  moves;  the  medium  within  each  cell  is 
represented  by  a  set  of  mass  points  (particles).  At  the  end  of  the  n-th  time 
cycle  the  mass  (equal  to  the  sum  of  the  masses  of  the  particles  located  in  that 
cell)  velocity,  pressure,  and  specific  internal  energy  are  associated  with  each 
cell.  To  obtain  the  corresponding  data  at  the  end  of  the  (n  +  1)  th  time  cycle 
one  makes  a  three  phase  calculation.  In  Phase  I  the  field  functions  are  changed 
neglecting  the  motion  of  the  medium.  In  Phase  II  the  mass  points  are  moved 
and  the  field  functions  then  recalculated  to  account  for  the  motion.  In  Phase 
III  various  functionals  are  computed  which  furnish  checks  on  the  accuracy  of 
the  calculations. 
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FORMULATION  OF  PROBLEM 


The  projectile  is  considered  to  be  a  circular  cylinder  of  length  L  and 
diameter  H.  An  axial  section  of  the  projectile  -  target  configuration  at  the 
instant  of  impact  (t  *  0)  is  shown  in  Figure  8;  it  is  superposed  by  the  fixed 
space  mesh  used  to  describe  the  subsequent  motion  of  the  configuration.  On 
this  plane  of  symmetry  the  cells  of  the  space  mesh  appear  as  rectangles  with 
sides  of  length  8r  =  h  by  Sz  =  k;  each  cell  is  actually  a  toroid  of  revolution. 
For  example,  the  projectile  is  initially  divided  into  /x  d  cells. 


/  =  L/k 


d  =  H/2h. 


(25) 


The  cells  are  labelled  with  index 
directions,  respectively: 


with  i  and  j  increasing  in  the  r  and  z 


cell 


(j-Dk< 
(i-Dh  < 


z<j  k 
r<i  h 


l  j 
*  i 


l,  2 . y 

1,  2,  •  •  • ,  / 3  . 


(26) 


Thus,  each  cell  includes  its  left  and  top  boundaries  but  does  not  include  its 
right  and  bottom  boundaries. 

The  projectile-target  material  is  represented  in  the  axial  plane  by  dis¬ 
crete  mass  points  called  "particles";  each  particle  is  actually  a  circle  about 
the  axis  of  symmetry.  Having  a  different  mass  for  different  mass  particles 
is  convenient  in  cylindrical  coordinates.  Each  particle  is  assigned  a  fixed 
mass  whose  value  is  proportional  to  the  radius  of  the  cell  within  which  it  lies 
originally;  i.  e. ,  at  t  =  0.  The  r  and  z  coordinates  of  each  particle  are  stored 
in  the  computing -machine  memory.  These  are  changed  in  time  in  accordance 
with  the  subsequent  motion  of  the  material  through  the  mesh  of  computational 
cells. 


If  there  are  originally  N  particles  in  each  cell,  then  the  mass  of  each  of 
the  N  particles  in  cell  is 


Z  ff  r ,  h  k  P 
1 _ o 

N 


2  ir(i  -  1/2)  h2  k  p 
_ o 

N 


a  (i  -  1/2)  h. 


(27) 


The  N  mass  particles  are  assumed  to  be  originally  randomly  distributed  within 
each  cell.  The  location  of  each  particle  at  the  end  of  the  n-th  time  cycle  is 
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stored  in  the  machine  memory  until  its  location  at  the  end  of  the  (n  +  l)th 
cycle  is  computed.  The  conservation  of  mass  is  therefore  automatic  in  the 
PIC  method. 

Such  dependent  variables  as  velocity,  pressure,  and  specific  internal 
energy  are  kept  in  the  machine  memory  by  cell.  For  example,  the  pressure 

of  cell  is  denoted  by  pj;  and  is  meant  to  represent  an  average  of  the  pres¬ 
sure  throughout  the  volume  of  the  medium  inside  the  cell  at  that  time.  The 
average  pressure  along  the  boundary  between  cells  and  is  denoted 

by  pj  while  the  average  pressure  between  cells  ^  ^  and  is  p1? 

Analogous  notations  are  used  for  the  other  cell-wise  quantities. 


DIFFERENCE  EQUATIONS 

The  velocity  components  and  the  specific  internal  energy  are  advanced  in 
time  according  to  a  finite  difference  approximation  of  the  equations  governing 
the  visco-plastic  model.  This  advance  occurs  during  Phase  I  of  the  calcula¬ 
tions  in  which  the  particles  are  not  moved;  thus,  the  transport  terms  are 
dropped  from  the  differential  equations.  The  tentative  new  cellwise  velocity 
components  become  (see  equations  A-2  and  A-3  in  Appendix  A): 


~j  „  „j  t  2  IMi-l/2>h  k  8t  \  pj  _pj  +*[  j  j  j 

i  i  kjj  |  i  +  1/2  i  -  1/2  h  [_  i  +  1/2  i+1  1 


-  ,  (u^  -  u^  ) 

i  -  1/2  '  l  i  -  l' 


+  ¥■ . 


j  4  it  (i- 1  / 2)h  k 


8 t  i  / u\  ^ 


m: 


j 


1/2 


(28) 


+ 


2  w  (i-l/2)h  8 1 


M' 


J 


a ! +  1/2 


1/2 


( 

I 
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-j  =  vj  +  Zrt(j.\IZ)h 
1  1  M 

i 


1/2  .  pj  -  1/2 


/2)h2  St  (  j  + 

i  t  ‘ 

1 

-f  -  W  *>]} 


(29) 


where 


-  a? 


j  j  +  1/2  j  -  1/2 

u.  V.  -  V. 


(i-l/2)h 


1  +  -L 


(30) 


2  ^ 

(D*)  =* 
i 


r  j  +  1/2  j  -  1/2  j  j 

i  i  l  +  1/2  l  -  1/2 

«  *  • 


+  2 


' 

Ui  +  1/2  “  Ui  -  1/2 

2 

X 

4  1 

i 

h 

T 

(i- l/2)h 

- 

- 

-I 


(31) 


(32) 


Pi  =  _  Pi  ‘  I  (diV  U) 


(33) 
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(34) 


Q?  ■  V3. 

1  i 


j  +  1/2  j  -  1/2  j  j 

u.  -  u:  v-  j.  i  o  ■  v*  i/o 

i  i  ^  i  +  1/2  x  -  1/2 


(r2/  «  (D2/  (Mp2. 

i  i 


(35) 


In  computing  the  tentative  new  cellwise  specific  internal  energy  both 


(u?#  v^)  and  (u1?,  v^)  are  assumed  available  in  the  machine  memory.  Then  with 

XI  11 


/ rr J 


u  =  |  (u  +  u) 


V  =  J  (v  +  v) 


(36) 


the  proper  form  is  (see  equation  A- 10  of  Appendix  A) 


-vj  j  +  2  ir  (i-  1  / 2)h  k  St 

{  *  (  y 


i  *  i 


l 


r  I  A? 


h  l'*i  +  1/2  '  Ai  -  1/2 


J 


Aj 
-  U. 
X 


r  *  (v  r 

rri+l/2  rr 


i+  1/2 


-(v__)J 


1  -  0j 
i~  1/2  J  1 

-  °!  [<v«>j + 1/2  -  <vj ' 1/2]^ 


i-  1/2 


(37) 


where  we  have  introduced  the  notations 


rj=,j+iihkl«  Aj^^j 


x 


(38) 


,i  =  i>  -  HJLMRi  Hi  pi  (dlvS,i 

1  1  M?  i 

X 


(39) 


19 


Ai  -  1/2 


=  -i  ,  (V  )j  +  $<V  )J+$  (V  )j 

2  1  -  1  rr  .  .  1  rr.  i-l  rz  .  , 

L  i-l  i  i-l 


j  -  1/2  1 


+ °f  <vrl>i] 

L  i  i  i 

(4 

+  0?  (V  )jl  . 
i  zz  , 

lJ 

The  viscous  stress  components  are  approximated  by  the  following  relations 
(see  equations  A- 11  of  Appendix  A): 


(V  ) 
rr 


j  =  aJ 

.  Mi  I  h  i  + 


AJ  . 

1/2  '  Ui  -  l/2) 


-4(divu)jj 


[(i-i%lh^  •  i (div^] 

<v„>;  • [|^tl/2-r  1/2i-|(di-,: 

<vr,»|  ■  ii  [i(fij  +  l/2-5|-  ‘/2>H 


1/2  Aj  -  1/2.  , 
-  u .  )  + 


±(£i  .0j  )1 

h  '  i  +  1/2  i  -  1/2'J 


wherein 


Aj  Aj  Aj  AJ  +  1/2  AJ-1/2 

“J  _  i  +  1/2  i  -  1/2  .  i  .  i  i 

(dlv  u). - - - +  __  + - - - 
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<D2)'  = 


A  j  +  1/2  Aj  -  1/2 

U  -  u 

i  i 


+  2 


,  +  1/2  '  -  1/2^ 

^  +  1/2  '  “i  -  I/2I  [  “i  1 
h  I  (i-  l/2)h  J 


Aj  +  1/2  Aj  -  1/2' 

V  -  V 

i  1 


-I  [«div"»i] 


S(I?,  p?) 

ir^l^y-rrtr-  ■ 

V  <  A  +  * 


(47) 


(48) 


PHASE  I  -  CALCULATIONS 

No  mass  particles  are  permitted  to  cross  the  left  boundary  of  the  mesh 
(axis  of  symmetry)  as  this  would  violate  the  assumption  of  rotational  symme¬ 
try.  No  such  restriction  applies  at  the  top,  bottom  and  right  boundaries  of 
the  mesh;  these  are  treated  as  "continuative  boundaries".  Accordingly,  the 
boundary  cells  along  these  three  sides  are  treated  as  interior,  being  bounded 
on  the  outside  by  cells  with  identically  the  same  properties  at  any  instant  as 
their  adjacent  interior  neighbors. 

Special  considerations  are  also  required  when  computing  in  a  cell  adjacent 
to  an  empty  cell  (if  the  cell  itself  is  empty  no  calculations  are  made).  The 
velocity  of  the  empty  cell  is  then  assumed  to  be  equal  to  that  of  the  cell  being 
computed;  the  pressure  and  the  viscosity  stresses  are  assumed  to  vanish  on 
the  boundary  of  an  empty  cell.  In  the  following  step-by-step  description  of 
the  computational  procedure  these  special  considerations  are  treated  in  de¬ 
tail  (a  programmer  may  work  directly  from  this  report). 

At  the  end  of  the  n  -  th  cycle,  the  mass,  specific  internal  energy,  velocity 
components,  and  the  pressure  for  every  cell  are  stored  in  the  memory.  Table 
IV  shows  these  together  with  the  quantities  which  replace  them  during  the 
sequence  of  computational  steps  required  for  Phase  I  of  the  (n  +  l)th  time 

cycle  calculations.  To  compute  a  given  step  for  cell  (A  the  results  of  the 
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previous  step  are  required  for  certain  cells  close  to 


(0- 


One  way  to  ensure 


this  is  to  require  that  each  of  the  computational  steps  is  made  for  every  cell 
in  the  mesh  before  the  subsequent  step  is  made  for  any  one  cell.  As  this  is 
the  simplest  method  to  de scribe ,  it  will  be  presented  here  even  though  it  is 
clearly  wasteful  of  computer  storage  space. 

Step  1.  1  The  velocity  components  at  the  end  of  the  n  -  th  time  cycle  are 
evaluated  at  the  cell  boundaries.  Ordinarily  one  uses  the  formulas 


uj 

1 


1/2 


1  l 


'J.i 


i  -  1/2 


(49) 


u1 


j  -  1  ,  j 

, , ,  u.  +  u: 

j  -  1/2  i _ i 


j  -  1/2 


i 


(50) 


If  cell  0  is  empty  at  this  instant.  *  0,  no  calculations  are  made  for 
this  cell  during  this  step  and  no  entries  are  made. 

If  cell  has  empty  neighbors  to  the  left  or  above  then  equations  (49) 
and  (50)  are  replaced: 


If  M?  ,  =  0  then  set 
i-l 


u. 


i- 1/2 


=  u? 


If  *  =  0  then  set 


j- 1/2 
u. 


v} 

i-1/2 


v?. 


l 


Empty  neighbors  below  or  to  the  right  cause  no  difficulty  in  this  step. 


If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  *  1  and  j  a  0,  1,  . . . , 
y  +  1,  then  replace  equation  (49)  by 


u 


1/2 


*i/a  *  i  (,vi  *  H'- 


The  first  relation  follows  from  symmetry  and  the  second  from  quadratic  ex¬ 
trapolation  with  (dv/dr)r_Q  =  0. 

9 

If  the  cell  is  adjacent  to  the  top  boundary,  j  *  1  and  1*0,  1,  ....  /3  +  1 , 
then  replace  equation  (SO)  by 


1/2 

v. 


No  special  provision  need  be  made  for  cells  adjacent  to  the  right  and 
bottom  boundaries. 

Step  1.  2  The  divergence  of  the  velocity  at  the  cell  center  is  computed. 
Ordinarily  one  uses  equation  (30). 

If  cell  is  empty  at  this  instant,  m|  =  0,  no  calculations  or  entries 
are  made. 


If  cell  ^ 


has  empty  neighbors  to  the  right  or  below  then  in  using  (30) 
one  must  use  special  provisions: 


If  *  0  then  set 


uJ  =  J 

Ui  +  1/2  V 


If  »  0  then  set 

l 


vi  +  l/2.vi. 

I  1 


* 
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No  special  provisions  are  required  for  cells  adjacent  to  the  axis  of  sym¬ 
metry  and  the  top  boundary. 

If  the  cell  is  adjacent  to  the  right  boundary,  i  *  0  and  j  »  0,  1 . 

Y  +  1,  then  use  (30),  but  set 


u 


0  +  1/2 


If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  i  *  0,  1,  .... 
0+1,  then  use  (30)  but  set 

y  +  1/2  =  y 

V.  *  V,  . 

1  1 

Step  1.  3  The  value  at  the  cell  center  of  the  distortional  strain-rate  in¬ 
variant  is  computed.  Ordinarily  one  uses  equation  (31). 

If  cell  is  empty  at  this  instant,  m|  *  0,  no  entry  is  made. 

If  the  cell  has  empty  neighbors  to  the  right  or  below,  then  in  using  (31) 
one  make 8  the  following  provisions: 

If  Mj+1  =  0,  set 


ui  +  1/2  =U^ 


'i  +  1/2  *  V 


If  =  0,  set 

i 


j  +  1/2  _  j 
u  *  u 

i  i 


j  +  1/2  j 
Vi  *  Vi* 


No  special  provisions  are  required  if  the  cell  is  adjacent  to  the  axis  of 
symmetry  or  the  top  boundary. 

If  the  cell  is  adjacent  to  the  right  boundary,  i  »  0  and  j  *  0,  1,  . . . , 

Y  +  1,  then  use  (31)  only  set 


u/9+  1/2  ’  ujs 
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If  the  cell  is  adjacent  to  the  bottom  boundary,  j  3  y  and  i  3  0,  1 . 

$  +  l,.then  use  (31)  only  set 


y  +  i/2  y 
i  i 


y  +  1/2 

i 


Step  1. 4  The  value  at  the  cell  center  of  the  strain-rate  dependent  vis¬ 
cosity  coefficient  is  computed  using  (32).  The  functions  if  (I.  p)  and  S  (I.  p) 
are  written  into  the  computer  code  as  subroutines. 


If  cell  ^  i 


is  empty  at  this  instant,  M?  3  0,  no  calculation  or  entry  is 


made.  Otherwise,  no  special  provisions  are  necessary. 

Step  1.  5  The  cell  center  values  of  the  quantities  P,  Q  and  are  com¬ 
puted  using  equations  (33),  (34)  and  (35). 


If  cell  is  empty,  3  0,  no  entries  are  made. 


If  the  cell  has  empty  neighbors  to  the  right  or  below,  then  in  computing 
Q  from  (34)  the  following  provisions  are  necessary: 

If  Mj+1  =  0,  set 

J  =  vi 

i  +  1/2  V 

If  Mj+1  =  0,  set 

l 

j  +  1/2  j 
u.  =  u;. 

i  i 

No  special  provisions  are  required  if  the  cell  is  adjacent  to  the  axis  of 
symmetry  or  the  top  boundary. 

If  the  cell  is  adjacent  to  the  right  boundary,  i  3  0  and  j  3  0,  1,  . . . , 
y  +  1,  then  the  calculation  of  P  and  cause  no  difficulty  but  in  using  (34)  to 
compute  Q  set 


v$  ♦  1/2  *  Vj9- 
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1, 


Similarly,  if  the  cell  is  adjacent  to  the  bottom  boundary,  j  3  y  and  i»0, 
. . .  ,  0  +1,  then  in  using  (34)  to  compute  Q  set 


u 


Y+  1/2 


i 


Step  1.  6  The  values  of  u/r,  P  and  p  at  the  left  cell  boundary  and  the 
value  ox  Q  at  the  top  cell  boundary  are  computed.  Ordinarily  one  uses  the 
formulas 


Mj  ,i  [-it.  t .  uL-1 

vv i-1/2  2  L(i  ■ 3/2)h  (i  ■  1/2)h-» 


(51) 


pj 

*i-l/2 


PJ  ,  +  PJ 
l-l  i 


(52) 


J  +  J 
j  -  1  +  **i 

^i-1/2"  2 


(53) 


QJ. 


j- 1/2  _  Qi  l+Qi 


(54) 


If  cell  is  empty  at  this  instant,  3  0,  no  entries  are  made. 

If  cell  has  empty  neighbors  to  the  left  or  above  then  equations  (51) 
through  (54)  are  adjusted: 


If  M?  ,  3  0  then  set 

l- 1 


Ui  1  *  Ui  *n  t0  comPute 


1/2  *  ° 


1/2 


j  _  j 

^i-1/2  ”  ^i* 
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If  *  -  0  then  set 
1 


Q 


j- 1/2 
i 


0. 


If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  -  1  and  j  *  0,  1,  . . .  , 
y  +  1,  then  replace  equations  (51),  (52)  and  (53)  by 


j 


=  0 


1/2 


If  the  cell  is  adjacent  to  the  top  boundary,  j  s  1  and  is0,  1,  . . .  ,  /3  +  1. 
then  replace  (54)  with 


1/2  1 
Q.  =Q.. 


No  special  provisions  are  required  if  the  cell  is  adjacent  to  the  right  and 
bottom  boundaries. 


j 

Step  1.  7  The  cell  center  value  of  the  tentative  new  radial  velocity,  u.,  is 
computed  using  equation  (28).  1 


If  cell  ^  i 
If  c.U  (j) 


is  empty  at  this  state,  m|  =  0,  no  calculation  or  entry  is  made. 


has  empty  neighbors  to  the  left,  to  the  right  or  below  then 
special  provisions  are  required  in  using  (28): 


If  M?  ,  =  0  then  set 

l- 1 
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If  .  »  0  then  set 
i+1 


PT 


i  +  1/2 


=  0 


J  -  J 

Mi  +  1/2  Mi 


u-  +  1  =  u- 
l  +  l  i 


(^)i+i/z’2’h  [pi7s+rrr7i] 


If  =  0  then  set 


Qj+1/2  =  0. 

i 


An  empty  neighbor  above  the  cell  causes  no  difficulty. 

If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  »  1  and  j  *  0,  1, 
y  +  1,  then  in  using  (28)  set 

If  the  cell  is  adjacent  to  the  right  boundary,  i  *  0  and  j  =  0,  1,  . 
y  +  1,  then  in  using  (28)  set 


pJ  s  pJ 

0  +  1/2  0 

j  _  j 

"0  +  1/2  "  M0 

us  + 1  ■  us 


0 


*1 

0  +  1/2  2h 


0-  1/2  0 


-1' 
+  1/2  J 


28 


If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  i  -  0,  1,  . . .  , 
/3  +  1,  then  in  using  (28)  set, 


Qr+in.ar 

I  i 


Step  1.  8  Compute  the  value  of  (rQ)  at  the  left  cell  boundary  and  the 
values  of  P  and  p  at  the  top  cell  boundaries.  Ordinarily  we  use  the  equations 


(rQ)Li/2  =  !  ‘/ao}] 


.  ,  ‘  1  +  pj 

-1/2  _  i _ L 


(55) 


(56) 


j  -  1  ,  j 
j-l/2  _  »*i _ !a 


(57) 


If  cell  is  empty,  =  0,  no  calculations  or  entries  are  made. 

If  cell  has  empty  neighbors  to  the  left  or  above  then  the  following 
provisions  are  made: 

If  M?  ,  a  0  then  set 
l-l 


<rQ);-i/2  =  o* 


If  M'?  *  *0  then  set 

i 


pj-1/23Q 

i 


j- 1/2  _  j 
Mi  "  Mi- 

Empty  neighbors  to  the  right  or  below  cause  no  trouble  during  this  step. 
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If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  »  1  and  j  s  0,  1 . 

y  +  1,  then  replace  equation  (55)  by 


<rQ) 


1/2 


0. 


Adjacent  to  the  top  boundary,  j  *  1  and  i  «  0,  1,  ...,/?  +  1,  replace 
equations  (56)  and  (57)  with 


p.172  =  p? 
i  i 


1/2  -  1 
“i  s  V 


No  special  provisions  are  required  for  cells  adjacent  to  the  right  and 
bottom  boundaries. 

Step  1 ,  9  The  cell  center  values  of  the  tentative  new  axial  velocity,  v?,  is 
computed  using  equation  (29).  1 

If  cell  is  empty,  3  0,  no  calculation  or  entry  is  made. 

If  cell  has  empty  neighbors  to  the  right,  above,  or  below,  then  spe¬ 
cial  provisions  are  required  in  using  (29): 

If  .  =0  then  set 

i+l 


(rQ  )i 


i+1/2 


0. 


If  *  *  0  then  set 


vj  '  1  =  vj. 
i  i 
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If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  1«0,  1.  . . . , 
fi  +  1,  then  in  using  (28)  set, 


y  +  1/2 


Step  1.  8  Compute  the  value  of  (rQ)  at  the  left  cell  boundary  and  the 
values  of  P  and  ft  at  the  top  cell  boundaries.  Ordinarily  we  use  the  equations 


(rQ):! 


-1/2 


_  h 
=  2 


(i  -  3/2)  Q?  ,  +  (i  -  1/2) 


Q? 

i 


(55) 


j-1/2 

Pi 


pj  -  1  +  pj 
i  i 


2 


1 


+ 


2 


(56) 


(57) 


M?  =  0.  no  calculations  or  entries  are  made. 

l 

neighbors  to  the  left  or  above  then  the  following 

provisions  are  made: 

If  M?  ,  *  0  then  set 
i-l 


If  cell  ^  i 
If  cell 


is  empty, 
has  empty 


If  M1?  *  =  0  then  set 
i 


J-1/2  . 


a  0 


j-1/2  _  j 
I-  =  Mi 


Empty  neighbors  to  the  right  or  below  cause  no  trouble  during  this  step. 
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If  *  0  then  eet 


p(+l/2.0 


1  1 

An  empty  neighbor  to  the  left 


of  cell 


cause  s  no  difficulty. 


If  the  cell  is  adjacent  to  the  top  boundary,  j  *  1  and  i  *  0,  1,  . . . ,  0  +  1, 
then  in  using  (29)  set 


0  1 

V.  a  V. 
1  1 


If  the  cell  is  adjacent  to  the  right  boundary,  i  =  0  and  j  *  0,  1,  . . . , 
y  +  1,  then  in  using  (29)  set 

(rQ>0  +  i/;  =  h0  Q@' 

If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  i*0,  1,  . . 
0+1,  then  in  using  (29)  set 


py+i/2  apy 

i  i 


y+1/2  y 

U.  =  U' 

1  1 


No  special  provision  is  required  for  the  cells  adjacent  to  the  axis  of  sym¬ 
metry. 
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Step  1.  10  Compute  the  cell  center  values  of  the  mean  velocity  components 
according  to  the  relations 


u *  -  {uj  +  u^) 
i  2  1  i  V 


If  cell  ^ 


(58) 


is  empty  at  this  stage ,  *  0,  no  calculations  or  entries  are 

made.  Otherwise  there  are  no  special  provisions  required. 

Step  1.  11  The  mean  velocity  components  are  evaluated  at  the  cell  bound¬ 
aries^  Ordinarily  one  uses  the  formulas 


Aj 

Ui  -  1/2 


AJ  .  Aj 
Ui  -  1  +  Ui 


$ 

i  -  1/2 


+Aj 
i  -  i  i 


(59) 


Aj  -  1  ,  Aj 

Aj  -  1/2  ^i _ IS 

Ui  2 


0j* 

05  -  1/2  _ Li 

i  2 


(60) 


If  cell 


is  empty,  *  0,  no  calculations  or  entries  are  made.  Other¬ 
wise,  the  same  special  provisions  that  are  detailed  in  Step  1.  1  are  required 
here;  only  with  u  and  v  replaced  by  Q  and  v. 

Step  1.  12  The  divergence  of  the  mean  velocity  at  the  cell  center  is  com¬ 
puted.  Ordinarily  one  uses  the  formula  (46). 

If  cell  is  empty,  3  0,  no  calculation  or  entry  is  made.  Otherwise, 

the  same  special  provisions  that  are  detailed  in  Step  1.  2  are  required  here; 
only  with  u  and  v  replaced  by  Q  and  0. 

Step  1.  13  The  cell  center  values  of  ¥  and  are  computed  using  for¬ 
mulas  (39)  and  (47). 


If  cell 


(0 


is  empty,  m|  3  0,  no  calculations  or  entries  are  made.  No 

other  special  provisions  are  necessary  for  V.  The  same  special  provisions 
that  are  detailed  in  Step  1.  3  for  the  computation  of  D2  are  required  here  in 
computing  82  except  that  u  and  v  are  replaced  by  Q  and  0. 


32 


Step  1.  14  The  cell  center  value  of  p  is  computed,  ordinarily  one  uses 
formula  (48). 

If  cell  is  empty  at  this  stage,  m|  3  0,  no  entry  is  made.  Otherwise, 
no  special  provisions  are  required. 

Step  1.  15  The  cell  center  values  of  Vrr,  Vrz  and  Vxz  are  computed, 
ordinarily  using  equations  (42),  (45)  and  (44)  respectively. 


If  cell  ||  j  is  empty, 

M?  3  0,  no  calculations  or  entries  are  made. 
i 

If  the  cell  has  empty  neighbors  to  the  right  or  below  then  the  following 
provisions  are  required: 

If  M? .  .  3  0,  set 
l+l 

Aj  _  Aj 

Ui+l/2  Ui 

when  using  (42) 

i+1/2  i 

when  using  (45). 

V  M?+1  3  0.  set 

l 

Aj+1/2  Aj 

Ui  =  Ui 

when  using  (45) 

0j+1/2 
i  i 

when  using  (44). 

No  special  provisions  are  required  if  the  cell  is  adjacent  to  the  axis  of 
symmetry  or  the  top  boundary. 

If  the  cell  is  adjacent  to  the  right  boundary,  i  3  fi  and  j  3  0,  1,  . . . , 
y  +  1,  then  set 

when  using  (42) 

when  using  (45). 
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If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  i  *  0,  1,  . . . , 
/3  +  1,  then  set 


Ay +1/2  a  y 

U  3  u 

i  i 


Ay +1/2  Ay 

V  =  V 


when  using  (45) 
when  using  (44). 


Step  1.  16  Compute  the  cell  center  value  of  Vqq  and  the  left  cell  bound¬ 
ary  value  of  A.  Ordinarily  one  uses  formulas  (43)  and  (40)  respectively. 

If  the  cell  is  empty,  =  0,  no  calculations  or  entries  are  made. 

If  cell  has  an  empty  neighbor  to  the  left,  i.  e. ,  ^  3  0,  then  re¬ 

place  (40)  by 


A3 


-1/2 


=  0. 


If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  *  1  and  j  3  0,  1 . 

y  +  1 ,  then  replace  (40)  by 


A  l/z  *  °- 


Here  we  have  taken  the  average  of  the  first  two  terms  and  the  average  of  the 
last  two  terms  of  (40)  to  vanish  independently,  since  the  radial  velocity  and 
the  shear  stress  both  vanish  at  the  axis  of  symmetry. 

No  other  special  provisions  are  required  for  \/^  and  none  at  all  are 
required  for  (V^  )j. 

Step  1.  17  Compute  the  cell  center  values  of  Vmax  and  r,  and  the  upper 
cell  boundary  value  ofO.  Ordinarily  one  uses  formula  (38)  for  I\  formula 
(41)  for  Q  and  the  following  formula  for  Vmax: 

{iw:i-  iw;i  •  i(v-)!i)- (6,> 

3  0,  no  calculations  or  entries  are  made. 


(Vm«)  [  ’ 


max 


If  the  cell  is  empty, 
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J-l 


If  cell  ^  has  an  empty  neighbor  above,  i.  e.  ,  Mj|~*  3  0,  then  replace 
(41)  by 


n 


j-  1/2 


0. 


If  cell  is  adjacent  to  the  top  boundary,  j  s  1  and  i  3  0,  1,  +  I, 

then  replace  (41)  by 


nl/2_Al  1  Al  .1 

u .  =  u.  (V  ).  +  v.  (V  ),  . 

l  l  rz  i  i  zz  i 


No  other  special  provisions  are  required  forQ^  and  none  at  all 
required  for  T  ?  and  V 


are 


max 


Step  1.  18  Compute  the  upper  cell  boundary  values  of  V  and  V  . 
Ordinarily  one  uses  the  formulas  rz  ZZ 


<Vrz> 


(V  )j'1  +  (V  )j 
.  ,  rz'  rz'. 

J-l/2  _  I _ i 


(62) 


.  W1  <vr-<v„>; 

<vzz>j  — —z — 1 

ZZ  Ct 

l 


(63) 


made. 


If  cell  i 


is  empty  at  this  stage,  3  0,  no  calculations  or  entries  are 


(V  )j'1/2  =  0 


has  an  empty  neighbor  above,  i.  e. ,  m|~  1  =  0,  then  set 


(V  )j‘1/2*o. 


zz  . 

i 


35 


If  the  cell  is  adjacent  to  the  top  boundary,  j  *  1  and  i*0,  1, 
then  replace  (62)  and  (63)  by 


(V  ) 
rz 


1/2 

i 


(V  ) 


rz 


1 

i 


<v»>1/2  ■  <v«>! 
i  i 


No  other  special  provisions  are  required. 

Step  1.  19  Compute  the  left  cell  boundary  values  of  V  and 
narily  one  uses  the  formulas 


(V_  ) 


(V  )J  +  (V  )J 
rri-l  rr  i 


rr 


i- 1/2 


(V  )J  +  (V  )J 
j  rZi-l  r*i 

<v  >J  3 — *-4 — L 

i-1/2  2 


If  the  cell  is  empty  at  this  stage,  M. 
made. 


0,  no  calculations  or 


If  the  cell  has  an  empty  neighbor  to  the  left,  i.  e. , 


0, 


<V  >J  =  0 
rr  i-1/2 


(V  )J  *  0. 


rz  . 


i-1/2 


.  •  • ,  fi  +  1, 


V  .  Ordi- 
rz 


(64) 


(65) 


entries  are 

then  set 
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If  the  cell  is  adjacent  to  the  axis  of  symmetry,  i  3  1  and  j  3  0,  1,  . . . , 
y  +  1,  then  replace  (64)  and  (65)  by 


(V,  ■  5  (9  v{  -  vb 

rr  ,/2  8  12 


1 


<Vrz» 


3  0. 


1/2 


No  other  special  provisions  are  required. 

Step  1.  20  The  tentative  new  specific  internal  energy  is  computed  using 
equation  (37). 

If  the  cell  is  empty,  3  0,  no  calculations  or  entries  are  made. 

If  cell  has  empty  neighbors  to  the  right  or  below  then  in  using  (37) 
one  must  make  special  provisions: 

If  M?  . ,  3  0  then  set 
l+l 


*1*1/2  =  ° 


(V  ) 
rr 


i+1/2 


(V  )J  3  0. 
rZ  i+1/2 


If  3  0  then  set 

l 

n  J+l/2  =  o 


(V  }j+1/z  =  0 
rz  . 

l 

(V  )j+1/2  =  0. 
zz  . 

i 
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•  •  •  9 


If  the  cell  is  adjacent  to  the  right  boundary,  i  *  fi  and  j  *  0,  1, 
y  +  1,  then  use  (37)  to  compute  I  but  set 


a  4 

A3+l/2  0 


(V  ) 
rr 


0+  1/2 


=  (V  )J 

rr/3 


(V  >J  3  (V  )J 
3+1/2  r  3 


If  the  cell  is  adjacent  to  the  bottom  boundary,  j  *  y  and  i  s  0,  1,  . . . , 
3+1,  then  in  using  (37)  set 


ny+1/2  »o  y 


(v  )y+1/2»(v  )Y 

rz  .  rz  . 

l  l 


(v)m/2.(v /. 

ZZ  .  ZZ  J 


Note,  if  a  negative  value  for  the  tentative  specific  internal  energy  is  com¬ 
puted  per  cell 

Step  1.  21  The  tentative  new  cellwise  velocity  components  are  recomputed 
using  the  formulas 


0 


by  the  above  procedure, 


••  j 
set  Ij 


0. 


Tij  »  2  Sj  -  u j 

i  i  i 


^  »  2  Oj  -  v[ 


U  c.U  (j) 


a 


is  empty,  M,  *  0,  no  calculations  or  entries  are  made. 


No  other  special  provisions  are  required. 
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PHASE  n  -  CALCULATIONS 


At  the  end  of  Phase  I  there  are  stored  in  machine -memory  fourteen  quan¬ 
tities  for  every  non-empty  cell  ^  .  Table  V  shows  these,  together  with 

the  quantities  which  replace  them  during  the  sequence  of  the  Phase  H  calcula¬ 
tions.  The  mass  particles  will  now  be  moved  and  this  effect  will  be  taken 
into  account  in  recomputing  the  velocity  and  specific  internal  energy  as  de¬ 
tailed  below. 

Step  2.  1  This  is  a  preparatory  step  in  which  the  tentative  energy  and 
momentum  are  totaled  over  all  the  particles  in  the  cell  at  the  end  of  Phase  I: 


1 


(66) 


i 


M. 


uJ 

i 


(67) 


z\  =  M?  v?. 

i  li 


(68) 


If  cell  j 


is  empty,  =  0,  no  calculations  are  performed  and  no  entries 


are  made. 


Step  2,  2  The  quantities  in  the  four  indicated  registers,  nos.  10,  11,  12 
and  13,  of  each  cell  are  replaced  by  zero.  This  preparatory  step  is  carried 
out  even  for  the  empty  cells. 


Step  2.  3  The  particles  are  moved;  new  values  for  the  total  cellwise 
mass,  energy,  and  components  of  momentum  are  computed.  To  perform 
this  step  one  goes  through  the  particles  (that  remain  in  the  system  at  this 
stage)  one-by-one,  computing  the  new  location  of  each  and  the  contributions 
of  each  to  the  new  values  of  the  cellwise  quantities.  The  contributions  are 
permitted  to  accumulate  to  the  final  total  values  of  M\  E',  R'  and  Z*. 


Thus,  suppose  the  particles  have  been  labeled  and  that  the  v  -  th  particle 
is  of  mass  and  its  coordinates  at  the  end  n  -  th  time  cycle  are  (r^  ,  ). 

The  substeps  associated  with  the  motion  of  the  v  -  th  particle  during  Step  2.  3 
are  shown  more  explicitly  in  the  schematic  of  Table  VI.  They  are  also  de¬ 
scribed  below.  After  the  calculations  associated  with  the  v  -  th  particle  have 
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been  completed  the  floating  storage  capacity  is  available  for  the  calculations 
associated  with  the  ( V  +  l)th  particle,  etc. 

(Substep  v  .  1)  The  indices  of  the  cell  in  which  the  v  -  th  particle  is 
located  at  the  end  of  the  n  -  th  time  cycle  are  computed  by  the  formulas 


[largest  integer  £  r^/h] 
[largest  integer  {  iy/k] 


(69) 


These  integral  values  are  then  stored  in  the  Floating  Storage. 

(Substep  v.  2)  The  velocity  components  of  the  u  -  th  particle  are  com¬ 
puted  by  the  "area  averaging"  formulas  presented  in  Appendix  B.  Their 
values,  u v  and  v„ ,  are  then  stored  in  Floating  Storage. 

(Substep  y  .  3)  The  new  position  coordinates  are  computed  by  the  formulas 


T'v  *  rw  +  u  v8t 

z'v  »  Zy  +  V  y$t. 


(70) 


The  values  ry  1  and  then  replace  the  values  r„  and  zy  in  the  u  -  th  Parti¬ 
cle  Storage. 

(Substepy  .4)  The  indices  of  the  cell  in  which  the  v  -  th  particle  is 
located  at  the  end  of  the  (n  +  l)th  time  cycle  are  computed  by  the  formulas 

[largest  integer  $  r'  /h  ] 

(71) 

[largest  integer  $c^/k]. 


These  integers  then  replace  rv  1  and  in  the  Floating  Storage. 
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(Sub step  v  .  5)  The  contribution  which  the  V  -  th  particle  makes  to  the 
total  mass,  energy,  and  momentum  in  the  cell  arc  now  computed  using 

the  formulas 


a  m 
1¥  l 


AE* 

i* 


. ,  m 

*'  v*} 


— :  E: 


M 


J  i 


m. 


m  u. 
1  v  1 


J*'  „ 


m 


AZ_  * - r 


m: 


rZ  *mu  vJ 
J  1  V  1 


(72) 

(73) 

(74) 

(75) 


If  both  0  <  i*<  0  +  1  and  0  <  j*<  Y  +  1  the  particle  still  remains  within  the 
mesh  and  these  contributions  are  then  added  to  the  quantities  already  present 
in  the  corresponding  memory  registers  (10,  11,  12  and  13  respectively)  of 

cell  .  If  either  i*  or  j*  violates  its  above  inequality,  the  particle  will 

have  left  the  mesh  and  no  entry  is  made  in  any  cellwise  storage. 

Step  2.  4  The  cellwise  density  and  the  velocity  components  at  the  end  of 
the  (n  +  l)th  time  cycle  are  computed: 


_ 1 _ 

2ir(i-l/2)h2 


k 


Also  insert  the  value  of  M? 

i 


into  the  first  register. 


(76) 


(77) 

(78) 
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If  cell  ^  j  is  empty,  i.  e. ,  *  0,  at  this  stage  no  calculations  or  en¬ 

tries  are  made.  No  other  prov'sions  are  required. 

Step  2.  5  The  cellwise  specific  internal  energy  at  the  end  of  (n  +  l)th 
time  cycle  is  computed: 


+  (vf  ) 


]• 


(79) 


j' 

If  the  cell  is  empty,  M*;  =  0,  no  computation  or  entry  is  made.  If  (79)  gives 

j' 

a  negative  value  set  Ij;  =  0. 


Step  2.  6  The  cellwise  value  for  the  tentative  new  pressure  is  computed 


tinder  the  temporary  assumption  that  no  material  fracture  occurs  in  cell 

(0 


during  the  (n  4-  l)th  time  cycle.  According  to  (19)  we  have  the  following: 


i  4»  4' 

pj  -Up3.,  l\) 


aMpf,  if) 


£  o 


U<{ 


<  o 


(80) 


If  cell 


is  empty, 


■  0,  no  entry 


is  made;  otherwise  no  special 


provisions  are  required.  The  functions  f(P,  I)  and  h(p,  I)  are  written  into 
the  computer  code  as  subroutines. 


Step  2.  7  The  cellwise  pressure  at  the  end  of  the  (n  4*  l)th  time  cycle  is 
computed.  According  to  hypotheses  (20)  and  (24)  we  have  the  following: 


pj'  =  0  if  pj  <  0 

l  i 


and  (V  )f  I  A- 

max  i  ■  «  J  * 

J  if  •  300]  ln*,/to 

( 

-  J 

*  p^  otherwise 


(81) 
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If  cell  j  is  empty,  M^  =  0,  no  entry  ie  made;  otherwiee  no  epecial  provi¬ 
sions  are  required.  The  fracture  criterion  (81)  is  also  incorporated  into  the 
computer  code  as  a  subroutine. 

At  the  end  of  Step  2.  7,  the  information  in  the  machine  memory  is  in  such 
a  form  that  the  computations  for  the  (n  +  2)th  time  cycle  may  begin  immedi¬ 
ately.  However,  first  making  some  subsidiary  calculations  is  helpful. 


PHASE  III  -  CALCULATIONS 

In  the  impact  problem  the  total  mass,  total  axial  momentum  and  the 
total  energy  of  the  projectile -target  system  are  rigorously  conserved: 

M  =  M’  =  Mass  of  System  (82) 

3  8 

Z  =  Z'  *4  H2  LP  V  (83) 

8  s  4  o  o 

E  =  E '  LP  v  Z.  (84) 

s  s  8  o  o 

These  quantities  and  the  total  radial  momentum  of  the  system,  R(,  will  be 
computed  and  monitored;  such  checks  serve  to  indicate  machine  or  coding 
errors.  * 

An  interesting  calculation  is  to  obtain  the  above  quantities  for  the  projec¬ 
tile  alone. 


M  »  M'  *4  H2  LP  ,  Z  ,  E  ,  R  , 
P  P  4  o  p  p  p 


and  the  location  of  the  projectile -target  interface: 


Maximum  value  of  the  axial  coordinates  of  all  those 
projectile  particles  that  lie  within  the  i-th  column  of 
cells  at  the  end  of  the  (n  +  l)th  time  cycle. 


(85) 


(86) 


*A  restart  procedure  will  also  be  written  into  the  program  so  that  a  machine 
error  at  any  stage  of  the  calculation  will  not  require  recomputing  from  the 
beginning  of  the  problem. 
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This  information  is  helpful  in  delineating  the  primary  flow  (during  which  the 
projectile  is  expended)  from  the  secondary  flow  (during  which  the  cratering 
mechanism  is  essentially  a  cavitation  process. ) 

a)  Total  System  Functionals 

To  compute  the  total  mass  of  the  projectile -target  system.  M#,  it  is 
only  necessary  to  cumulate  the  masses  of  the  particles  within  the  mesh  and 
the  masses  of  those  that  have  left  the  mesh  by  crossing  one  of  the  three  con- 
tinuative  boundaries.  Each  of  the  departing  particles  will  also  carry  its  con¬ 
tributions  to  Es.  Rs.  Z#  outside  the  mesh.  In  addition,  however,  energy  and 
momentum  are  lost  across  the  continuative  boundaries  by  diffusion  once  the 
disturbance  produced  by  the  impact  is  propagated  to  the  mesh  boundary. 
Hence,  to  monitor  Eg,  Rg,  and  throughout  the  computation  it  is  necessary 
to  provide  for  loss  across  the  boundaries  both  by  particle  motion  and  by  dif¬ 
fusion: 


M'=(M').  +(M') 
s  s  in  s  out 


E'  *  (E1 )  +  (E' )  .  where  (EM  *  (E')di“*  +  (E')***1, 


s  in  s  out 


s  out  s  out 


s  out 


R'  »  (R'K  +  (R1)  .  where  (RM  _  *  (R')di**  +  (R')*** 


s  in  '  s  out 


s  out  s  out 


s' out 


Z'  .  (Z'l.  +  <Z'>  .  wh.re  {Z-) _ »  +  (ZM^f 


s  in  '  s'out 


s  out  s  out 


s  out 


(87) 

(88) 

(89) 

(90) 


The  schematic  in  Table  VII  depicts  the  processes  by  which  the  above 
quantities  are  cumulated  from  the  contributions  carried  by  the  particles  and  by 
diffusion  from  cells  just  inside  continuative  boundaries.  The  contributions 
carried  by  the  v-th  particle,  whether  it  remains  within  or  departs  from  the 
mesh  during  the  (n  +  1)  th  time  cycle,  are  shown  explicitly.  The  particle 
contributions  are  actually  cumulated  during  Step  2.3.  The  losses  by  diffusion 

of  E'  R(  Z(  across  the  boundary  of  cell  ^  ^  ,  with  i»/8orj*lorj»y, 

are  denoted  by 


A  {  (E^) 


,  \8iff* 


out 


Aj  (R.)di«- 

i  s  out 


<z;C 


(91) 
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The  total  diffusive  losses  are  given  by 


(  ) 


diff. 

out 


I 

Cells  at 
Cont.  Bdry. 


4i< 


diff. 
out  * 


(92) 


where  the  contributions  (91)  are  computed  by  the  formulas  presented  in  Ap¬ 
pendix  C.  The  contributions  to  the  diffusive  losses  of  momentum  and  energy 
are  actually  cumulated  during  Steps  1.  6  and  1.  16.  respectively,  as  indicated 
in  Table  VII. 

These  calculations  are  therefore  not  made  subsequent  to  Phase  n  but  are 
cumulated  during  both  Phases  I  and  II  of  the  time  cycle.  At  the  start  of  each 
time  cycle  the  register  in  which  the  contributions  to  (  )^n  are  to  be  cumulated 
will  be  zero,  but  each  register  in  which  the  contributions  to  (  )out  are  to  be 
cumulated  may  be  non-zero,  since  it  includes  the  losses  by  parti.de  motion 
and  by  diffusion  during  the  previous  n  time  cycles. 

At  the  end  of  the  (n  +  1)  th  sweep  through  the  cells  and  the  particles,  that 
part  of  each  functional  which  remains  within  the  mesh,  (  )jn,  and  that  part 
which  has  left  the  mesh,  (  )out,  are  entered  individually.  In  the  last  substep 
shown  in  Table  VII  these  corresponding  subtotals  are  combined  according  to 
(87) . (90). 


b)  Delineation  of  Projectile 

The  functionals  (85)  and  (86)  are  of  interest  only  during  the  first  stages 
of  the  cratering  process.  Therefore,  no  provision  for  diffusive  losses  will 
be  made  and  these  functionals  will  lose  their  significance  in  the  later  time 
cycles. 

The  corresponding  cumulation  process  for  the  projectile  material  is  de¬ 
picted  in  the  schematic  of  Table  VIII.  In  this  case  one  rejects  all  particles 
of  the  target  material  and  then  cumulates  the  contributions  individually  ac¬ 
cording  to  whether  the  particle  remains  within  the  mesh  or  is  outside.  Finally, 
the  corresponding  contributions  are  combined  to  give 


M'-(M').  +(M') 
p  p  in  p  out 

(93) 

E«  *  (E' )  +(E')  , 

p  p  in  p  out 

(94) 
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R'  a  (R‘ ).  +  (R') 
p  p  in  p  out 

(95) 

Z'  a  (Z')s  +  (Z') 
p  p  in  p  out 

(96) 

Moreover,  the  projectile-target  interface,  defined  by  (86),  is  computed  dur¬ 
ing  this  sweep  in  the  manner  depicted  in  Table  VIII;  T).  denotes  the  maximum 
new  axial  coordinate  of  all  the  first  (v  -  1)  particles  tkat  lie  in  the  i  -  th  col¬ 
umn  of  cells. 


c)  Axial  Shock  Wave 


Along  the  axis  of  symmetry  the  solution  is  identical  to  that  of  the  impact 
between  two  semi-infinite  bodies  (under  corresponding  conditions)  until  rare¬ 
faction  waves  arrive  from  the  edge  of  the  projectile-target  interface.  The 
comparison  of  the  location,  amplitude,  and  velocity  of  the  shock  front  with 
the  one -dimensional  solutions  of  Reference  6  gives  a  further  check  on  the 
accuracy  of  the  calculations  during  the  early  stages  after  impact. 


CHANGE  OF  NET  SIZE 

As  time  goes  on  the  size  of  the  crater  increases  and  the  stress  wave 
propagates  further  into  the  target.  More  target  material  must  then  be  covered 
by  the  calculation  mesh  than  is  necessary  at  earlier  times.  As  the  dimensions 
of  the  disturbances  increase,  however,  sufficient  resolution  may  be  obtained 
by  using  a  larger  net  size,  in  both  time  and  space,  than  was  permissible  dur¬ 
ing  the  initial  stages  of  the  process.  It  will  therefore  become  advantageous 
to  repartition  the  system  during  the  course  of  a  computational  run. 


a)  Calculation  of  Time  Step 

The  length  of  the  time  step  Jt  may  be  varied;  a  smaller  value  is  required 
during  the  early  stages  of  the  cratering  mechanism  than  at  a  subsequent  time 
when  the  particle  velocity  has  decreased.  In  the  present  program  it  is  com¬ 
puted  for  each  step  from  the  data  of  the  previous  time  step.  The  idea  is  to 
specify  an  integer  as  the  maximum  number  of  particles  that  will  move  a  dis¬ 
tance  so  large  that  at  the  end  of  the  time  interval  the  particle  is  not  in  an 
adjacent  (or  identical)  cell.  The  specified  integer  represents  the  maximum 
number  of  times  that  condition  (B-l)  of  Appendix  B  may  be  violated. 
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b)  Repartition  of  Space  Mesh* 


In  rezoning  the  space  mesh  it  is  convenient  to  choose  the  original  number 
of  particles  per  cell,  N,  to  be  a  perfect  square.  Since  the  optimum  value  of 
N  is  apparently  between  3  and  9t  it  is  natural  to  set 

N  a  4 

in  the  present  part  of  the  formulation  since  an  increase  in  the  courseness  of 
the  mesh  by  a  factor  of  nine  may  well  be  too  abrupt.  An  additional  advantage 
of  the  choice  N  ■  4  is  that  for  the  same  machine -memory  capacity  it  permits 
smaller  values  of  h  and  k  which  in  turn  decrease  the  effect  of  the  artificial 
dissipative  terms  inherent  to  the  method  (Reference  33). 

Assume  that  in  the  existing  space  mesh  there  are  an  even  number  of  col¬ 
umns  and  rows  and  that  originally  there  were  four  mass  particles  in  each  of 
the  cells,  i.  e.  ,  /3  and  Y  are  even  and  N  =  4.  Also  assume  that  sufficient 
machine  storage  space  has  been  allotted  to  accommodate  4(3  x  Y  particles. 
Then  an  enlarged  space  mesh  may  be  introduced  in  which  the  linear  dimen¬ 
sions  of  the  cells  are  doubled,  i.  e. ,  four  of  the  original  cells  are  combined 
to  form  a  single  cell  of  the  new  mesh  (see  Figure  9).  The  area  covered  by 
the  space  mesh  may  thereby  be  increased  fourfold  without  increasing  the 
number  of  cells  in  the  mesh.  The  time  increment  may  also  be  increased  in 
the  new  mesh: 


h  -*  2h  k  -*  2k.  (97) 

In  the  75%  of  the  mesh  area  which  has  been  newly  introduced,  each  cell 
is  assigned  four  mass  particles  randomly  distributed  within  the  cell.  The 
mass  of  each  of  these  particles  is  given  by  (27)  with  substitutions  (97)  in¬ 
serted,  i.  e. ,  the  mass  of  each  is  eight  times  those  of  a  particle  initially 
(t  *  0)  in  one  of  the  original  mesh  cells  with  the  same  radial  index  i. 

In  the  remaining  25%,  which  is  being  repartitioned,  the  masses  in  the 
enlarged  cells  are  computed  in  the  following  manner.  Each  quadrant  (sub¬ 
cell)  of  an  enlarged  cell  is  precisely  one  of  the  original  cells.  The  mass 
particles  in  each  sub-cell  are  coalesced  to  form  a  single  particle  provided 
the  sub- cell  is  not  empty;  if  the  sub- cell  is  empty  before  repartitioning  the 


♦The  properties  of  the  projectile  alone  can  no  longer  be  isolated  after  repar¬ 
titioning  as  some  of  the  coalesced  masses  will  contain  both  projectile  and 
target  material.  The  calculations  depicted  in  Table  VIII  are  thereby  mean¬ 
ingless  after  repartitioning. 

tSee  the  discussion  following  equation  (108). 


— — 


corresponding  quadrant  of  the  enlarged  cell  remains  empty.  The  position  of 
the  coalesced  mass  particle  is  taken  to  be  the  mass  center  of  the  particles  in 
the  corresponding  sub-cell  prior  to  coalescence: 

r«  *  Ma  Y  r-  <98) 

sub- cell  a 

z«  =  jl  52  ***  zv  (M«  =  S  >• 

Cl 

sub-cell  a  sub-cell  a 


The  values  of  Ma  are  available  from  the  last  cycle  of  calculations  made  prior 
to  repartitioning. 

The  velocity  of  the  mass  particles  in  the  newly  introduced  75%  of  the  net 
area  may  be  assigned;  it  is  zero  in  the  present  problem.  The  velocity  of  each 
of  the  coalesced  masses  in  the  original  25%  of  the  area  may  be  chosen  in  such 
a  way  as  to  conserve  energy  and  momentum.  Conservation  of  axial  and  radial 
momentum  requires  that  the  corresponding  velocity  components  of  the  en¬ 
larged  cell  be  given  by 


v 


M 


4 

X  Ma  va 
a  a  1 


(100) 


u 


M 


M«  u« 


4 

(M  *  X 

a  = 


Ma). 


(101) 


The  ua  and  va  denote  the  sub-cell  velocity  components  and  are  available  from 
the  last  cycle  of  calculations  made  prior  to  repartitioning. 
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To  conserve  total  energy  the  specific  internal  energy  of  the  enlarged  cell 
must  be  computed  from  the  relation 


.  2.  ,  i 

+  v  >  +  M 


4 


1 

a  =  1 


+  7<u 


2 

a 


,  2.  .  1 

tv  )  +  M 


I 


(102) 


where  the  Ea  denote  the  total  energy  in  each  of  the  sub-cells  composing  the 
enlarged  cell;  these  values  are  available  from  the  last  cycle  of  computations 
made  prior  to  repartitioning.  The  density  of  material  in  the  enlarged  cell  is 
computed  from 


P 


M  > 
h  k  r 


(103) 


Here  h  x  k  are  the  dimensions  of  the  original  net  and  r  denotes  the  radial 
distance  to  the  center  of  the  enlarged  cell.  The  pressure  in  the  enlarged 
cell  is  then  computed  according  to  the  equation  (according  to  relations  15 
and  20): 

p  *  i(p,  I)  if  C  £  0 
*  h(p,I)  if  C  <0 
9  0  otherwise. 

From  (100)  and  (101)  one  may  calculate  the  decrease  in  the  total  kinetic  en¬ 
ergy  within  the  enlarged  cell: 


•nd-j-  £  (V  )  < 
4  «  .  1  mix« 


cr 


(104) 


AK.  E.  =  y  M  (u2  +  V2)  -  Y  i  Ma(ua+Va) 
c  a  =  1 


_  i 

—  X  M  M  ' 

2M  a  a 
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Hence, 

AK.  E.  ^  0. 


49 


(105) 


There  is  a  loss  in  kinetic  energy,  AK.  E.  <  0,  except  for  the  case  in  which 
the  velocities  of  the  four  sub-cells  are  identical.  However,  equation  (102) 
insures  that  the  total  internal  energy  in  the  enlarged  cell  is  increased  by  the 
same  magnitude. 

Actually  the  total  radial  momentum  of  the  system  is  not  conserved.  This 
suggests  an  alternate  procedure  in  which  (91)  and  (92)  are  replaced  by  the 
conditions  that  the  kinetic  energy  and  internal  energy  are  individually  con¬ 
served: 


u 


M«  <u«  +  v«> 


(101*) 


M«  la¬ 


st  =  l 


(102*) 


In  this  case  the  total  radial  momentum  in  the  enlarged  cell  is  increased  by 


AR.  M. 
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=  Mu  -  2 

a  a  1 


M« 


s  a,  a  ?  l  M«  M«  [  (u«  -  u«  >2  +  <va  -  V  )2 

4  ,)l/2  4 

+  <  i  Ma  ua>  1  -  «?  i  M«  u«  • 


Hence, 

A  R.  M.  Z  0. 


(103') 


Equality  holds  only  if  the  velocity  of  each  of  the  four  sub-cells  is  the  same. 

Which  of  these  procedures  gives  the  more  accurate  results  will  have  to 
be  decided  by  actual  computational  runs.  In  either  base  the  computational 
method  detailed  in  the  preceding  paragraphs  may  be  carried  out  for  subse¬ 
quent  time  cycles  using  the  enlarged  mesh.  It  is  only  required  that  substi¬ 
tutions  (97)  be  made  in  the  equations  governing  each  step  of  the  calculations. 
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Mass,  energy,  or  momentum  that  has  previously  departed  from  the 
mesh,  ..either  by  particle  motion  or  by  diffusion  across  the  mesh  boundaries, 
cannot  be  recovered  by  the  repartition  process  described  above.  Thus,  the 
quantities  that  should  be  rigorously  conserved  during  repartitioning  are 


(M  ).. 
s  in 


(E  >.  . 

s  in 


<Z  hn 
s  in 


rather  than 


M  ,  E  ,  Z  . 
s  s  s 

If  the  repartitioning  is  delayed  too  long  the  loss  of  M,  E,  R,  Z  across 
the  continuative  mesh  boundaries  will  introduce  large  errors.  To  avoid  this 
the  computer  program  may  be  written  to  automatically  repartition  whenever 
the  cumulative  loss  of  any  one  of  the  four  quantities  exceeds  a  specified  small 
fraction  of  the  total  quantity  remaining  inside  the  mesh. 
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COMPUTER  PROGRAM  (PICWICK) 


The  computational  procedure  described  in  the  preceding  sections  taxes 
both  the  memory  capacity  and  the  speed  of  most  computers.  The  particular 
programming  logic  selected  depends  upon  the  trade-off  between  computation 
time  and  storage,  whichever  is  most  advantageous  for  the  particular  computer 
used.  The  following  is  a  brief  description  of  the  considerations  that  have  led 
to  PICWICK,  a  computer  program  for  the  IBM  7090  that  uses  only  internal 
storage.  No  account  has  been  taken  of  the  possibility  of  using  a  high-speed 
magnetic  disc  for  intermediate  storage  since  none  is  presently  available  at 
either  the  General  Electric  Company  or  at  the  Eglin  Air  Force  Base.  The 
use  of  magnetic  tape  would  be  prohibitively  slow. 


STORAGE  FOR  CELLS 


To  compute  a  given  step,  in  the  Phase  I  and  Phase  II  calculations,  the 

results  of  the  previous  step  for  certain  cells  near  to  must  be  available. 

From  a  bookkeeping  standpoint  this  is  most  simply  accomplished  by  requiring 
that  each  computational  step  be  made  for  every  cell  in  the  mesh  before  the 
subsequent  step  is  made  for  any  one  cell.  This  straightforward  procedure, 
depicted  in  Tables  IV  and  V,  would  require  a  total  of  14  0Y  storage  locations 
for  the  cellwise  quantities.  This  number  would  be  required  for  the  Phase  I 
■calculations  and  would  also  provide  for  the  Phase  II  storage  by  "overlays". 
This  simple  procedure,  however,  is  wasteful  of  machine  storage  space. 


The  required  cellwise  storage  can  be  greatly  reduced  by  carefully  se¬ 
quencing  and  by  making  use  of  movable  storage.  If  computing  time  was  of 
minor  importance,  the  storage  requirements  could  be  reduced  to  5/8  y  loca¬ 
tions.  A  method  has  actually  been  wbrked  out  which  accomplishes  this  by 
trading  computing  time  for  storage.  This  breaks  down  as  follows:  5  locations 
per  cell  are  needed  in  Phase  I  to  store  M,  u,  v,  I  and  Vmax,  with  T*  "read 
out"  in  Phase  I  when  called  for.  Modifications  in  the  criterion  for  fracture 
could  reduce  this  to  4  locations  by  eliminating  the  need  to  store  Vmax.  In 
Phase  II,  5  locations  are  needed  to  store  M,  u,  v,  I  and  p.  The  method 
would  probably  not  be  economical  with  the  IBM  7090,  but  probably  would  be 
with  the  IBM  7030  ("Stretch"). 


To  save  time  on  the  7090  in  the  Phase  II  calculations  an  additional  three 
locations  per  cell  are  required  for  AE,  A  R  and  A  Z.  Hence  the  total  is  eight 
per  cell.  Since  the  storage  can  be  shared  with  Phase  I  the  total  cellwise 
storage  required  for  the  PICWICK  program  is  8/8  y . 
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STORAGE  FOR  PARTICLES 


For  each  of  the  N 0y  particles  one  location  is  needed  to  store  u,  and  one 
to  store  v.  Since  there  are  only  N/S  distinct  masses  all  the  masses  can  be 
stored  in  N/3  locations  by  proper  grouping.  Therefore  2N/9  y  +  N/3  *  2N/3y 
locations  are  apparently  required  for  the  particles.  Unfortunately,  however, 
books  must  be  kept  on  the  particles  that  have  left  the  mesh.  To  do  so  re¬ 
quires  N/3  locations,  which  can  be  "shared"  with  the  N  locations  required  for 
the  particle  masses. 

Therefore,  3N/3y  locations  are  needed  for  the  particle  storage  if  particles 
leave  the  mesh  and  approximately  2N/3Y  are  needed  if  they  do  not.  The  for¬ 
mer  is  applicable  to  the  present  problem. 

TOTAL  PERMANENT  STORAGE 

The  storage  required  for  the  program  itself  is  quite  large.  For  compar¬ 
ative  purposes,  however,  the  necessary  total  permanent  storage  is  approxi¬ 
mated  by  the  sum  of  the  cell  and  particle  storage  requirements  for  a  P  x  Y 
cell  mesh  with  originally  N  particles  per  cell.  If  the  simple  procedure  of 
keeping  books  is  followed  the  storage  requirement  would  then  be  approximated 
by 

Sj  ■  (14  +  3N)0Y  (106) 


The  more  sophisticated  bookkeeping  procedures  used  in  the  PICWICK  code 
reduce  the  total  permanent  storage  required  to 

S2  «  (8  +  3N)/9y.  (107) 


In  order  to  obtain  some  estimate  of  the  optimum  number  of  particles  per 
computational  cell,  in  a  memory  limited  computer  such  as  the  IBM  7090,  the 
heuristic  argument  of  Harlow  (Reference  35)  may  be  applied.  Accordingly, 
one  optimises 


Q^/SY^N  and 


0Sy)3/2N 


(108) 
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with  S  held  constant.  The  first  alternative  gives  equal  weight  to  the  Eulerian 
and  Langrangian  systems;  the  second  choice  gives  a  greater  weight  to  the 
Eulerian  system.  Applied  to  Sj.  the  calculation  yields  14/3  for  the  first 
alternative  and  28/3  for  the  second.  Applied  to  the  results  are  8/3  and  8 
respectively. 

The  optimum  number  of  particles  per  cell  is  apparently,  therefore,  from 
N  s  3  to  N  *  9.  For  most  of  our  calculations  the  choice  of  N  *  4  is  used  as 
we  have  seen  this  to  be  especially  convenient  in  repartitioning  the  space 
mesh  -  a  desirable  procedure  if  the  cratering  process  is  to  be  studied  through¬ 
out  its  formation.  If,  indeed,  four  particles  per  cell  are  used,  the  IBM  7090 
internal  memory  (32,  000  registers)  would  be  exhausted  by  a  mesh  of  f$  x  y  < 
1230,  about  30  by  40  cells,  if  the  simple  bookkeeping  procedure  leading  to 
(106)  is  used.  On  the  other  hand,  if  the  program  storage  requirement  is 
approximated  by  (107)  a  mesh  of  /3  x  Y  <1600,  about  35  by  45  cells  would 
become  feasible.  The  latter  applies  to  the  PICWICK  program.  This  is  very 
important  because  the  nearly  optimum  resolution  thus  obtained  is  apparently 
just  sufficient  to  make  useful  calculations  possible  without  resorting  to  super 
computers  such  as  Stretch. 
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CONCLUSIONS 


A  "debugged"  computer  program  of  the  visco-plastic  governing  equations 
has  been  developed  for  axisymmetric  hypervelocity  impact  between  a  pro¬ 
jectile  and  target  of  the  same  material.  The  computational  scheme  requires 
the  equation  of  state,  the  flow- resistance  coefficient,  and  the  fracture  criteria 
for  dynamic  tri -axial  stress  conditions.  Since  none  of  these  phenomenological  re  - 
lations  have  been  firmly  established  by  experiments  under  the  severe  conditions 
operative  during  hyper  velocity  impact,  the  proposed  relations  are  considered  tenta¬ 
tive,  and  have  been  written  into  the  program  as  sub-routines  so  th  at  they  may  be 
readily  changed.  The  uncertainties  involved  are  less  important  than  the  ability  of 
the  analysis  to  include  and  compare  various  proposed  phenomenological  relations. 

In  the  special  case  that  the  flow  resistance  coefficient  is  set  equal  to  zero 
and  the  effect  of  material  strength  is  ignored,  the  governing  equations  reduce 
to  the  perfect  fluid  equations.  In  the  first  computer  production  runs,  to  be 
made  with  PICWICK  at  Eglin  AFB,  this  case  will  be  treated.  This  will 
provide  a  check  with  the  calculations  presented  by  Bjork  (Reference  38). 

A  realistic  theory  must,  however,  include  the  effects  of  flow- resistance 
and  fracture  resistance.  In  the  first  runs  in  which  these  effects  are  included 
i)  and  S  will  be  assigned  constant  values.  This  will  permit  comparison  of  the 
results  with  those  previously  obtained  for  the  case  of  normal  impact  between 
two  semi -infinite  iron  bodies  (Reference  6).  During  the  early  stages  after  impact 
the  axisymmetric  solution  along  the  axis  of  symmetry  should  be  in  agreement 
with  the  one -dimensional  solution  obtained  using  the  same  values  of  *70,  SQ, 
and  impact  velocity.  Subsequent  runs  will  be  made  for  copper  and  aluminum 
using  the  estimates  (13)  for  q0  and  handbook  values  for  S0. 

The  early  production  runs,  in  which  the  dependence  of  ft  on  the  thermo¬ 
dynamic  state  is  included,  will  investigate  the  effect  of  the  temperature  while 
neglecting  the  effect  of  the  pressure.  There  are  two  reasons  for  this  pro¬ 
cedure.  More  experimental  data  are  available  concerning  the  effect  of  temper¬ 
ature  on  the  flow- resistance  of  materials  than  the  effect  of  pressure.  Further, 
hypervelocity  impact  experiments  in  which  the  temperature  of  the  target  ma¬ 
terial  was  widely  varied  have  shown  that  drastically  different  crater  dimen¬ 
sions  result;  crater  volumes  of  lead,  cadmium,  zinc,  and  copper  were  ob¬ 
served  to  increase  by  a  factor  of  two  to  five  when  the  target  temperature  was 
near  the  metal's  melting  point  (Reference  39).  Indeed,  the  inability  of  the 
perfect  fluid  model  to  account  for  these  results  was  one  of  the  factors  that 
led  us  to  develop  the  visco-plastic  model  of  hypervelocity  impact.  In  these 
computer  runs  the  forms  (8)  and  (14)  will  be  used  in  which  the  parameters 
are  approximated  respectively  from  plastic  wave  propagation  experiments 
and  static  yield  tests  carried  out  at  elevated  temperatures. 
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APPENDIX  A 


Let  p,  p,  u  =  (u,  0,  v)  and  I  denote  the  density,  pressure,  velocity,  and 
specific  internal  energy  respectively.  Then  the  axisymmetric  formulation  of 
the  vis co -plastic  equations  may  be  written  in  the  form  (References  2  and  6): 


M  +  “lf  +  vlf  + » ,uv  a  ■  ° 

*  (if +  u  It + vls) =  £  (p  +  2  p-  It) 

+  2  P  <!•  P‘  D)  £  (7)  +  If 

*  (lr  +  u  It +  v  lt)*lr  (p  +  2  p-  D>ti) 

^(^■+u^“  +  v  +  P  div  u  *  M(I»  p,  D)  D2 


P  -  F  (p,  I) 


where 


P  =  "  P  *  "5  M  (I,  p,  D)  div  u 
Q  =  Md,  P.  D)  (||  +  |j) 

(•  .  P  (I,  p,  D)  ■  V  0,  P)  +  ,  S(I,p) 

iT7T 7 

r2  =  D2  p2 


(A-l) 


(A-2) 


(A-3) 


(A-4) 


(A-5) 


>  (A-b) 


D2  . 


(st  +  It)  2  +  2  [(It)  2+(t)  2+  (£)  *] 


•|  (dlvu)  2 

divu  =  t*  t  & 
or  r  d* 
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The  introduction  of  the  parameter  e  >  0  into  the  definition  of  the  vis¬ 
cosity  coefficient,  p  (I,  p,  D)  greatly  simplifies  the  problem  since  it  elim¬ 
inates  the  moving  surface  of  separation  between  the  rigid  and  fluid  regions  of 
the  medium.  Thus  is  not  necessarily  >  r  This  not  only  simplifies 
the  calculations,  but  is  also  more  realistic  because  the  shearing  stresses 
now  depend  upon  the  strain- rate  in  a  continuous  manner  (See  Figure  3). 


For  improved  stability  of  the  numerical  calculations  it  is  advantageous 
(Reference  35)  to  write  the  right  side  of  the  energy  equation,  (A-3),  in  an 
alternate  form.  First  note  that# 


M  (I,  p,  D)  D2  =  \  V  *.  D, 


(A-7) 


where  D  is  the  deformation  tensor  and  V  is  the  viscous -stress  tensor.  Upon 
using  the  identity 


div  (V  .  u)  =  (div  V)  .  u  +  j  V  :  D  (A-8) 

it  fpllows  that  (A -4)  may  be  written  in  the  form 

p(li +  °  t?+  v  tr)  +  p  div3  ■ div(*  • 31  • (div^  • a-  (a-9) 


If  the  right  side  is  expanded  some  of  the  terms  cancel  and  one  obtains 


#In  Reference  6,  $  and  fi  were  denoted  by  and  respectively.  In 
Reference  37  the  latter  quantity  is  denoted  by  2D. 


61 


where  the  components  of  the  viscous  stress  tensor  are  given  by 


V__  *  M  (D _ -  |  div  Q)  *  M  (2  -  |  div  tl) 


rr 


rr  3 


%  =  M  (D00  -  |  div  n)  ■  (2  j  -  |  div  u) 

VZZ  *  »  <Dzz  -  4  dlv  «>  ■  *  <2 

V  =  n  D  =  M{|^  + 
rz  rz  oz  or 


(A 


The  normal  components  of  7  are  the  deviatoric  stresses  and  Vrr  +  V, 


+  V„  =  0. 


-ID 
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APPENDIX  B 


Consider  a  particle  in  cell  .  During  the  (n  +  1)  th  time  cycle  it  may 
remain  in  that  cell  or  it  may  move  into  one  of  the  nearby  cells.  The  particle 
would  necessarily  remain  within  its  current  cell  or  move  into  one  of  the  eight 
nearest  neighbor  cells  if  were  restricted  to  satisfy  the  inequalities. 


St 


u 

max 


<  h 


St 


V 

max 


<  k. 


(B-l) 


The  area  averaging  procedure  described  below  does  not  require  this  restriction, 
and  it  has  been  found  advantageous  to  relax  the  condition  to  permit  (B-l)  to 
be  violated  by  several  particles  at  each  step.  This  permits  fit  to  be  chosen 
larger  and  greatly  reduces  the  computation  time  with  insignificant  loss  of 
accuracy. 

Imagine  a  rectangle  of  cell  size  located  about  each  particle,  the  particle 
being  at  the  center  (see  sketch).  The  effective  velocity  for  moving  the  particle 
is  taken  as  the  weighted  average  of  the  four  cellwise  velocities  of  the  cell  that 
the  superposed  rectangle  overlaps.  The  weightings  are  proportional  to  the 
overlap  areas. 


Consider,  for  example,  the  v  -th  particle  and  suppose  it  is  currently 
located  in  cell  The  appropriate  values  of  i  and  j  as  well  as  the  current 

coordinates  of  the  particle,  (r^  ,  z  y),  are  available  at  the  end  of  Substep  y .  2 
(see  Table  VI).  Now  the  cells  which  the  superposed  rectangle  overlaps  vary 
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with  the  quadrant  of  cell  (  J  )  within  which  the  particle  lie*.  Each  pos¬ 
sibility  and  its  corresponding  tentative  new  cellwise  velocity,  8y  ■  (ay,  vy), 
are  described  in  the  following: 

First  Quadrant 

(i-l/2)h-  ry  <  ih 
(j-l)k^  ay  <(j-l/2)k 

%•'-  I  [<j-l/2)k  -  «v](  [r„  -  (i-l/2)hj  *  [<i+  l/2)h  -  r  J  ‘‘) 


[ 


]  ([« + »«»> -  %]  ’{ *  [  %  -  +i) 


+  (j-3/2)k 

4  \  U  rj  *  b  -  j 

(B-2) 

If  one  of  the  overlapped  cells  is  empty  special  provisions  are  required  in 
using  (B-2): 


IfM^'j  =  0  then  set 


..j-1  -j 
ui+l  i 


If  m|”1  *  0  then  set 


-r  -  *t 


If  *  0  then  set 


am  ■  °Ji 


4 
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If  cell  is  adjacent  to  the  right  boundary  of  the  mesh,  i 

j  =  0,  1 ,  . . . ,  y+1,  then  use  (B  -2)  only  set 


~j-l  _ -j-1 

u3+l  '  u/5 


-j-1  -j-1 

3+i  3 


-j  - 

u3+i  ”  a3 


3 +1  3 


If  cell  ^0  is  adjacent  to  the  top  boundary  of  the  mesh,  j  =  1 
0+1,  then  use  (B-2)  only  set 


-0  -1 
Ui+1  "  ui+l 


-0  _1 
Vi+1  '  vi+l 


-0  _  _1 
Ui  ‘ai 


-0  -1 
Vi  “  Vi 


Second  Quadrant 


(i  -  l)h^  rv<(i  -  l/2)h 
(j  -  l)k^  zy<(j  -  l/2)k 


\=~j[(j-1/2)k-Zv]  ([rv-^-3/2>h]  ai'l+[(i’ 
<]  ^(i-l/2)h-ryJ  _j+^ 


■  1 /2)h  - r 


+  -  (j  -  3/2)k  |  [|(i-l/2)h-rv| 


rv-(i-3/2)h 


1 


If  one  of  the  overlapped  cells  is  empty  special  provisions  are 
in  using  (B -3): 


If  M1?  *  =0  then  set 
i 


uj_1  =uj 
i  i 


d-i  s  J 

vi  vi* 


0  and 


and  i  =  0, 


required 
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If  m|"|  b  o  then  set 


■  5t 

i-l  l 

t|  j  ■  0  then  set 

-j  _  - j 
ui-l  “  ui 

V*  -  V^. 

i-l  i 

If  cell  is  adjacent  to  the  top  boundary  of  the  mesh,  j  *  1  and  i  »  0, 

. . ,  0+1,  then  in  using  (B-3)  set 

n°  .  s1 

-0  .1 
tr  M 

®i  Ui 

Vi  Vi 

-0  .1 

Ui-1  *  Ui-1 

-0  -1 

V  »  v  * 

i-l  i-l 

If  cell  y  )  is  adjacent  to  the  axis  of  symmetry,  i  ■  1  and  j  ■  0,  1,  , . , , 
-  1,  then  musing  (B-3)  set 

_j-l  _j-l 

Ur  *  .«** 

0  1 

0  1 

^ 

0  1 
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Third  Quadrant 


(i-  l)h^  ry  <  (i-  l/2)h 
(j  ■  l/2)k  “  Z|/  <  jk 


If  one  of  the  overlapped  cells  is  empty  special  provisions  are  required 
in  using  (B-4): 


If  M 


j 

1-1 


0  then  set 


~J  _  ~J 
Ui-1  '  Ui 


If  M^+j  =  0  then  set 
l-l 


~j+l  _~j 
Ui-1  "  Ui 


v^1^. 

i-1  i 


If  M^+1  =  0  then  set 
i 


~j+l  ~j 
Ui  =Ui 
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If  cell  ^3^  U  adjacent  to  the  axis  of  symmetry,  i  ■  1  and  j  *  0,  1,  , . . , 
Y  +  1,  then  in  using  (B-4)  set 


uo ■ -  “l 


j+1  j+1 

uo  =  -ui 


V*  S  v3 
0  V1 


=  v3*1 
0  V1 


If  cell  ^3  ^  is  adjacent  to  the  bottom  boundary  of  the  mesh,  j  ■  y  and 
!■  0,  1,  .,,1  0+  1,  then  in  using  (B-4)  set 


~Y  +1 

~y 

~y+l 

~y 

Vl 

•  Vi 

Vl 

B  V 

Vi 

-wy.J 

Ui 

~y 

8  If 

-.y+l 

-  vY 

Ui 

vi 

Vi  * 

Fourth  Quadrant 


!(l-i/2)h<  vy  <ih\ 

(j  -  l/2)k<  *y  <jk  ) 

u  -i.  ||(j+l/2)k-*,J  ^  -(i-l/2)hj5j+1+  £(i+l/2)h-^J  S3^ 

+  [■  -  (j-  l/2)kj  ^  ju+i/2)h-ryJ  53+1+  £r  -  (i  -  W2)hJ  53* 


(B-5) 


If  one  of  the  overlapped  cells  is  empty  special  provisions  are  required  in 
using  (B-5): 

If  M3+  0  then  set 


u3  *  u3 

Vl  i 


v3  -  v3 

i+1  V 
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If  M^=  0  then  set 


~j+l  _  ~j 

ui  "  ui 


~j+l  _  ~j 

V  -  V 

i  i‘ 


If  =  0  then  set 

l+l 


~J+ 1  ~j 
u. , .  =  u. 
i+l  l 


If  cell  Q  is 


~j+ 1  ~j 

V.  =  V.. 
1+1  1 


adjacent  to  the  bottom  boundary  of  the  mesh,  , 
i  =  0,  1,  ....  3+1,  then  in  using  (B-5)  set 


~Y+1  ~y 
u.  =  u. 
i  i 


~Y  +  1 
v.  =  v. 
i  i 


~y+l  ~y 
ui+l  =  ui+l 


— y+1  ~y 

V  •  V 

i+l  i+l 


If  cell  is  adjacent  to  the  right  boundary  of  the  mesh,  i 

j  =  0,  1  ....  y  +1,  then  in  using  (B -5)  set 


uj  =  uj 

3+1  3 

~j  ~j 

v3+i  "  v3 

~  j+ 1  ~  j+ 1 

u3+l  =  u3 

l> 

II 

l  > 

j  =y  and 


=  3  and 
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APPENDIX  C 


In  the  computational  procedure  the  normal  and  shearing  streeeee  along 
each  continuative  boundary  are  asaigned  to  be  juet  equal  to  their  correspond* 
ing  values  in  the  adjacent  interior  cell.  Surface  forces  are  thereby  assigned 
to  the  mesh  once  the  disturbances  produced  by  the  impact  are  propagated  to 
the  mesh  boundary.  The  momentum  and  energy  produced  by  the  surface 
loading  must  be  accounted  for  if  the  conservation  relations  (83)  and  (84)  are 
to  remain  valid  throughout  the  computations.  For  each  cell  along  a  con¬ 
tinuative  boundary  it  is  necessary  to  make  the  following  additions  to  the 
energy  and  momentum  cumulatione: 


Right  Mesh  Boundary  (excluding  corner  cells): 


i 

=  fa  and 

j  =  2,  3,  . .  •  • 

y  -  i 

A  i(E'  )diff* 
/?  s  out 

n 

c  > 

[p  j  -  (V  )j 
L0  rrp. 

]  -0p«vr.>pj[2'f* 

Aj  (R‘  )dif 

fa  8  OUt 

■-p  j, 

^2*  fthk  gt  j 

A*  (z*  )diff 

fa  s  out 

■-Q  0 

|^2ir  0hk  gt  j 

Bottom  Mesh  Boundary  (excluding  right  corner  cell): 


j  s  y  and  i  *  1,  2,  . . . ,  fa -l 

Af<E,.C =  |  n  [-! ]  -  sr  <v„>r  |  [2»<*->'2>  ** »«] 

ar«R,.c^rK‘^> 1,2  »*] 

Af(z'.C=-pi  [2'<‘-‘/2)  h2»*] 
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Top  Mesh  Boundary  (excluding  right  corner  cell): 

j  B  1  and  i  *  1.  2,  . . . , 0  -1 

[A -«v„>11]+ai1«v„),lj[2,(i-,/2»  h2»‘] 

A‘i  <R’.>"“  ■  Ql  »2»<] 

4!<z’.C  -pi  [*»<i-i/*>h2»t] 

Lower  Right  Mesh  Corner,  i  s  fj  and  j  *  Y  : 

^,e’.c  I2**1**] 

♦jy  [»s  •  <;  (vr.>i}[2-o-i/2>  »2»'] 

^<R'.C  *  -Pp  [z»  Phkjtl  -Q^2»0-l/2)h28tJ 

^(Z'.,o«  '  'QP  [2'  Phk*t]  -PP  [2’<0-1/2>  h2#*] 

Upper  Right  Mesh  Corner,  i  *0  and  j  *  1: 

*0  ,E'.C  fy  [-r  <vrr)^]  >0  <v„)is}  [2^1*  «.] 

+{u0(v„>e-C  lt  [»*r  <v„)‘a] }  [*•»-»/*)  »2**] 
4<«'.C  ■  -pjj[j»**at]  +Q/3 [2«-0-x/2)  h2«t] 
a3(z'<1^  ■  *Q a  [**<&*•* ] +  pjb  [2»(/3-i/2i  h2»t] 
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formation  is  Reference  13 


TABLE  II 


Equation  of  state  parameters,  occurring  in  relation  (15). 
expressed  in  the  gm-cm*  c  system  of  units. 


Aluminum 

Iron 

al 

1. 1867 

7.  78 

a2 

0. 76300 

31.  18 

b 

3. 4448 

9.591 

o 

bi 

1. 5451 

15.676 

b2 

0.  96430 

4.634 

c 

0.43382 

0. 3984 

o 

C1 

0. 54873 

0. 5306 

*o 

1. 5000 

9.  00 

2.702 

7.86 

£  j  =  =  0  for  t  <  t* 

£  =  £  =  0.  045  for  t  >  t* 

X  Cm 

Here  t*  is  the  time  after  which  a  major  signal  has  passed 
over  the  material. 
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TABLE  HI 


Values  of  the  material  parameters  required  in  the  strength  -  time  relation, 
equation  (24)  used  to  predict  fracture  (References  22  and  24). 


Xx  10'4 

lrt-7 
e  x  10 

t 

o 

Metal 

(degrees) 

(degrees/mb) 

(£sec) 

Aluminum  (annealed  at 

2.  55 

5. 1 

2.5  x  10‘5 

550°C) 

Aluminum  (annealed  at 

2.5 

3 

2.5  x  10-5 

420°C ) 

Zinc 

1. 25 

0.8 

6. 6  x  10 

Silver 

3.  2 

0.9 

3.  55  x  10'8* 

Nickel 

4.  35 

0.3 

1. 55  x  ID'7" 

Platinum 

6.  25 

3.85 

7.9  x  IQ’7 

*  Obtained  from  single  0  vs  log  t  curve  at  18°C 
**  Estimated  from  formula 

'oSh<kV‘‘ 

where  h  is  Planck's  constant,  k  Boltzmann's  constant  and  ^  D  is  the  Debye 
temperature  (Reference  25). 
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STEP  1 


TABLE  IV 
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STEP  1 


** 

d 

& 

fc 

0, 

fc 

& 

0 

u 

W 

W 

W 

W 

W 

U 

H 

H 

H 

H 

H 

H 

w 

Vi 

(0 

Vi 

w 

V) 
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Sequence  of  cellwise  storage  changes  for  cell 
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STEP  2.  7 


a;  A  A  A  A  A 

4) 

q  ^  ^4  oO 

3  i  i  i  *  £ 

w 


79 


TABLE  VII 


Schematic  of  the  cumulation  processes  for  the  phase  III  calculations 
associated  with  the  projectile -target  system. 


M 

E 

R 

z 

s 

* 

1 

I 

0 

0 

0 

0 

— 

— 

— 

— 

Add 

Add 

Add 

Add 

0 

0 

0 

0 

At  End  of  n-th  Tim*  Cycle 
At  Start  of  (n  +  l)th  Tim*  Cycl* 


(M  )  „ 

•  out 

(E  )  „ 

•  out 

‘Vout 

<Z*>out 

It 

II 

II 

II 

Diffusive  Loss  el 

[Stop  1 . 6) 

Add 

Add 

Add 
i  Ml 

*K, 

Add 

Contlnuatlv* 

0 

0 

Boundary  of... 
Math:  C*U  ({) 
with  1  *0or 

1  *  1  or  j  *  y 

[St*p  1.16 

Add 

M 

Add 

Add 

0 

0 

0 
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projectile 


Figure  1.  Illustration  of  projectile -target  configuration  just  before  impact. 


-5  MSEC  b.  15.  6  MSEC 


e.  54.  6  MSEC  f.  89.  8  m SEC 


Figure  2.  Illustration  of  violent  ejection  of  material  from  target  during  the 
formation  of  a  crater.  In  (a)  the  projectile-target  system  is  de¬ 
picted  just  prior  to  impact;  the  subsequent  ejection  is  depicted  in 
the  sequence  (b)  to  (e).  (Sketched  from  photographs  obtained  at 
BRL) 
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Figure  3.  Schematic  representing  the  dependence  of  the  shearing  stress 
on  the  rate  of  deformation:  (a)  Effect  of  c ,  (b)  Effect  of  in¬ 

ternal  energy  (or  temperature)  and  pressure. 
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Figure  4.  Viscosity  factor  depicted  as  a  function  of  the  inverse  of  the 

absolute  temperature,  ijis  given  in  the  gm-cm-/isec  system 
of  units  (from  References  11  and  12).  The  curves  are  broken 
for  temperatures  in  the  solid  phase  range. 
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Figure  5.  Schematic  of  equation  of  state: 

(a)  Function  equal  to  pressure  in  compressive  region 

(b)  Extrapolation  into  tensile  region. 


Figure  6.  Strength- time -temperature  behavior  of  platinum  (References  22 
and  24):  (a)  Dependence  of  time  to  fracture  on  the  absolute  tem¬ 
perature  (b)  Dependence  of  X  -  a to  on  o. 
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Figure  7.  Strength- time -temperature  behavior  of  aluminum  (References  22 
and  24)  (a)  Dependence  of  time  to  fracture  on  the  absolute  tem¬ 
perature  (b)  Dependence  of  X  -  u>a  on  o . 
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TOP  CONTI  NUATIVE  BOUNDARY 


mm 


(i)  (!)  (!) 


d  =  NO.  OF  CELL  COLUMNS  IN  PROJECTILE 
1--  NO.  OF  CELL  ROWS  IN  PROJECTILE 

NO. OF  COLUMNS  IN  MESH 
y  -  NO. OF  ROWS  IN  MESH 


(D  (2)  (3) 


BOTTOM  CONTINUATIVE  BOUNDARY 


Figure  8.  Rectangular  mesh  superimposed  on  the  projectile-target  con¬ 
figuration  at  instant  of  impact  (t  =  0). 
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Figure  9.  Schematic  representation  of  the  repartitioning  in  which  the  mesh 
area  is  increased  four-fold. 


9 


APGC  -  TDR  -62-74 


INITIAL  DISTRIBUTION 


2 

Wpns  Sys  Eval  Gp 

4 

Bu  of  Nav  Wpns  (RM) 

1 

Hq  USAF  (AFTAC) 

1 

Nav  Rsch  Lab 

1 

Hq  USAF  (AFCIN-3K2) 

2 

Naval  Ordnance  Test  Station 

2 

Hq  USAF (AFRDC) 

2 

Naval  Ordnance  Test  Station 

1 

Hq  USAF  (AFRAE-E) 

(Tech  Lib) 

1 

Hq  USAF  (AFTST-EL/CS) 

2 

Naval  Ordnance  Lab  (Tech  Lib) 

1 

Hq  USAF  (AFRST-PM/ME) 

2 

Naval  Wpns  Lab  (Tech  Lib) 

1 

AFSC  (SCRWA) 

1 

USAF  Proj  RAND 

1 

AFSC  (SCTA) 

1 

NASA 

2 

BSD  (AFSC) 

4 

NASA  (Langley  Rsch  Ctr) 

1 

ASD  (ASAPRL) 

4 

NASA  (Ames  Rsch  Ctr) 

1 

ASD  (ASRMDS) 

1 

Marshall  Space  Flight  Ctr 

3 

ASD  (ASAD- Lib) 

(M-AERO-TS) 

1 

ASD  (ASRNGW) 

1 

Marshall  Space  Flight  Ctr 

2 

SSD 

(Advanced  Rsch  Proj  Lab) 

1 

OOAMA  (OOYD) 

1 

Advanced  Rsch  Projects  Agcy 

2 

AFSWC  (Tech  Info  Div) 

1 

Defense  Rsch  and  Engineering 

1 

AFCRL  ( CRQST-2) 

(Tech  Lib) 

2 

AFOSR 

2 

Defense  Rsch  and  Engineering 

3 

AFOSR  (SRHP) 

1 

Armour  Rsch  Foundation 

1 

AFOSR  (Univ  of  Utah) 

(Def  Rsch) 

1 

TAC  (TPL-RQD-M) 

3 

Lewis  Rsch  Ctr 

2 

ARO  (Scientific  Info  Br) 

1 

Univ  of  Chicago  (Institute  for 

2 

OAR  (RROSA) 

Air  Wpns) 

1 

White  Sands  Msl  Range 

2 

The  Franklin  Institute  of  the 

(ORDBS-OM-W) 

State  of  Penn  (Tech  Lib) 

3 

AMC  (MCR) 

1 

Calif  Institute  of  Technology 

2 

Picatinny  Arsenal 

(Jet  Propulsion  Lab) 

(ORDBB-TK) 

2 

Johns  Hopkins  University 

4 

Aberdeen  Proving  Ground 

(Applied  R  sch  Lab) 

1 

Redstone  Scientific  Info  Ctr 

30 

ASTIA  (TIPCR) 

1 

Frankford  Arsenal  (Lib) 

APGC 

1 

Frankford  (Pitman-Dunn  Lab) 

2  ASQR 

2 

Springfield  Armory  (R&D  Div) 

2  PGAPI 

2 

Watervliet  Arsenal  (ORDBR-R) 

3  PGEH 

2 

Watertown  Arsenal 

2  PGWR 

1 

Rock  Island  Arsenal 

20  PGWRT 

1 

Engr  Rsch  Sc  Dev  Lab 

2  PGW 

(Tech  Support  Br) 

1  USCONARC  (Spc  Wpns  Div) 

2  U,  S,  Army  R sc h  Office 
4  Bu  of  Nav  Wpns  (R-1Z) 


Page  91  of  91  pages 


I  a  £  x  i  s 

!l.y*S3*  •* 

Lit  m; 


au‘  g 

SH  5 


-~"*-ag  >: 


H  jJ  t 

O  5  y  * 

phi 

u.  Ik  •  *•  ii 

.•  6  h  i 3 

m* 


!|S| 

mQ  >  o  • 
U  »  «  Z 

-  . «  *  2 

t  *•  ft.  ,  - 

£  *j  Z  4  _ 

u  ,  u2  j 

*s  *  s  *.  5 

-Saifs 

mis 

o  .  3  "  « 

**rr 

5  io  i* 

4***o- 


y  C  *3  "5  C  ^  *  £  f* 

is  « 3  *  SE 

£§gijsl»-i!siasi2e|-g 

1  i *  | *  Ii sf  f  ■*  •  •  5 !  s  b  “  i 

i!  H  li!  pill! !° 

>.-8^|x;laS«s|6f«Jft4 

"saJ 


IlitiHli 


•  t  i  3  s ;  s  a  .s  4 1  S  |  s  s  a  *  - « 

IU  UH  Stf  f  ? :  5 1 1|! 

i il: *11 54 Ills  Is  s  il 


l  nu  i 

\ui\m 

si  a  >  > 


i  l 
? 


0*3 

1SH 

!'is‘l= 

nu 

SSSs 

§  1 
58S2- 
t*wr:3 
iiali 

1? } 

m 

iuti 


iLthll 


&  5  *  *  c  S  B  4  2  s  i  3  i  3  8  £  e  5 

.}|sS|a Ijtlslfi 

ili  hi !  j  ii 


koiJ!SSH:  J3uS{ 

IliliSSfiil:! 


dillsisilfUmm 


_  > 

,il, 

ml  I 


imm 

i;B.fcouo 


'  a  g  >  > 


•  **'  e 

3  QSi 

i  th  i 

lllsslKl! 

a  g  >  > 


P  s 

I  Iff 

Uh 

IgjtS 
hS  0  I  S 
3  2  2  •- 

•Oh!3 

ii** 

"  >  o.  . 

s&’s 

ss§i 

nil 

Sin- 

..In! 

tzi’t 

2?*IS- 
18fi|l 
t»S|  | 

?i&ia 

1 13"  j 

£  3  <  t  5 

sS84* 

<  *  *  *  o- 


s  •  -  4  , 
si .  s I •= 

•  *8  3  8  2  ' 

>  s.  i  o  -•  2  ' 
S  0  5  o  >-  S  , 

iljfii 

"Sis  .  j 

ii?Hi 

*  llifi 

i;  • :3  * 

2 %s  . Sa- 

a  I  IS Sf 
3  1  s  -s  j??  . 

lilslf 

:is£sl; 

111  s'  |1 

»  w  e  «  S  ;  ; 

»■ « s  *?  * 

y  .  ;  f  •  »< 

"  !  s  «  h  »  ’ 
c  *  5  O  ^  > 

!  Ill 


r;  2  ■»  •  *  a  g  „ 

•  Ss'3  iJS 

°  -  f  2  2  I  3  •  S  8  C  f  1 


its silica! sit 


HU 


ifcH-SfSg,!! 
l85el®<8S'If*S 
1 1  x  i  I  P-  2  jj  ?  S  S  1 J 

«  J  U  N  y  6  a.  C  a  I  w9m 


a  “  «* 

jHi 


!#5* 

lit-. 

?sj!i 

u  .  u£  J 

i?!if 

«ih 

s  .8"  * 

Anti 

<***  o- 


dlteSHW!0 

Ha  ii  ii 


m 


;  *  u  rf  3 

if'S:  3 

**  I.  y  9  ri  3 


'  m  "  *  -  a  a  >*• 


8  tit 

5sU 

O  B  U  OC 


f  §  -i  I  3  |  I  si.  £ 

t’sll ;l<  s  s ;  *!-  .? ;  as 
i  S  S 1  i  5  S  •  1 1 8 1  *  a  { ? 


o*  a 

3  gat* 

iIJm 

fceji  *■  * 
•  6  hi 3 

|8SJ 

*;* 

■*?£ : 

5PX^ 

1 3s  f 
58  >Ji 

■  .  jj  "*  s 
I  *76-* 

I  N  “  i#  . 

UTUl* 

18S|I 

VI* -z 

£lTr 

±  s.o  2*  * 

<  «  i*  m  o- 


II  f  Is  2  2 1  sfi 

Huflbiiivag  iftfi  D.I  U3S 


tj  si?. 

4  s»<3“;  i 

*  SfSl  * 

is  ilii 


lEiff 
35 S  la 
c  E 1  s  X 

*0fc'3D 

\ii* 

Ull 

fp] 

58>L 

•  .  3  *•  C 
i*  ♦  E  .  * 

HI 

Hi!| 

55  Ik 


pi  1  fimiiio 

111!-.  !  Ifslil-Jl 


Hit  i"{ IHitll 

u  8  M  lit!  !  Ini'!  ! 


s  1 1 4  s 
IHVI 


J  •  u  *.  w*2o  "» 

J SaS‘4 ! 

ifsaSI ill* 

:;uoui!  >3 


“ a  3 


|  e  -  : 

7,  2,  •«*  *  » 


H  *X  »• 

<  ;  e  l 
9  2  S  « 
O  £  Z  K 

nh 

3g?«3 

z  zt  o  a  • 

224  N 

:§gl3 


-23* 

l  11 . 

u  Sjli 

»  N  "  J  . 

u’*u£  » 
v  g« 

IBfil  a 

f  5 

f &2s  a 

r.s-i 

£  |  <  t  a 

S  &8  §•  i 

<  et  u.  jc  , 


j  i  i  -1  £  .  2 

8 1  a  s  t  i c  ?  | 
•3l«S‘  .o! 


H-«S!  *•  «s  ‘  :  - 


Siilu 

.  o  w  •o  • 

s  •  i°  •  • 

Hi  4*  I 

5  A  1. 2  3  | 

hil»  3 


0t.S56!.‘l 

.■^ttOuuS  aa 

J{-j  O.^J 

l|l|l  Hi 

3 il *  is Ilf 


Hill 

3Un-iVjs!»l 

a  i  .5  , 1  5  s  <  3  S  i  s  « 
!•  S;!!-1"’” 

3|§l35jE!i?*S 
s  stlflf s I  it 
i  5  5  1 1  ?g  !  jt  s  F.?g 

\j  v  a  c  a.  >  w 


2  2 1 

?*,s .  I 
ifli  I 


Iliallll-I 


■aa  >:> 


§H! 

Jl# 

131m 

hfc*  ‘  K 

lisT 


WiWwKi^ 

tijHlljs  }  |1  iiO 

It I i||f iil| 

f  in  fiisill \m 


JJBtJ 

JSoS  .• 

1?  ?l 


isi 

f  I  ll 

ijtSI} 


i!  ii*  f  |i?!l 

! rli Is  S  ll 


